R/FolderAnalysis_functions.R
In CNVPanelizer: Reliable CNV detection in targeted sequencing applications

Documented in Background BedToGenomicRanges BootList CNVPanelizer CNVPanelizerFromReadCounts CNVPanelizerFromReadCountsHELPER CollectColumnFromAllReportTables IndexMultipleBams NormalizeCounts PlotBootstrapDistributions ReadCountsFromBam ReadXLSXToList ReportTables SelectReferenceSetByInterquartileRange SelectReferenceSetByKmeans SelectReferenceSetByPercentil SelectReferenceSetFromReadCounts StatusHeatmap WriteListToXLSX

###############################################################################
# functions to read the read counts
###############################################################################

BedToGenomicRanges <- function(panelBedFilepath,
                               ampliconColumn,
                               split = "_",
                               doReduce = TRUE,
                               rangeExtend = 0,
                               dropChromossomes = NA,
                               skip = 1) {
  #load the bed file
  segments <- read.table(panelBedFilepath,
                         sep = "\t",
                         as.is = TRUE,
                         skip = skip)
  gr <- GRanges(segments[, 1], IRanges(segments[, 2], segments[, 3]))

  #make sure that all regions are more than one read apart
  GenomicRanges::start(gr) <- GenomicRanges::start(gr) - rangeExtend
  GenomicRanges::end(gr) <- GenomicRanges::end(gr) + rangeExtend

  if(doReduce) {
    reducedGr <- GenomicRanges::reduce(gr, ignore.strand = TRUE)

    #get the genes for each reduced range
    hitsAnnotation <- GenomicRanges::findOverlaps(reducedGr, gr,
                                                  select = "first",
                                                  ignore.strand=TRUE)
    #get the amplicon names from the segments
    amplicons = segments[hitsAnnotation, ampliconColumn]
    #replace the original gr with the reduced gr
    gr <- reducedGr
  } else{
    amplicons = segments[ , ampliconColumn]
  }

  #get the genes from segments
  splitted = strsplit(amplicons, split = split)
  # Required because of package checking complaints
  i <- NULL
  genes = foreach(i = seq_along(splitted)) %do% {
    splitted[[i]][1]
  }

  genes = unlist(genes)
  # missing the dollar symbol before text because of groovy script
  elementMetadata(gr)["geneNames"] = genes
  elementMetadata(gr)["ampliconNames"] = paste(seqnames(gr),
                                               start(gr),
                                               end(gr),
                                               sep = "_",
                                               amplicons)

  if (!is.na(dropChromossomes)) {
    gr <- dropSeqlevels(gr, dropChromossomes)
  }
  return(gr)
}

#load the files you want to analyze
ReadCountsFromBam <- function(bamFilenames,
                              sampleNames,
                              gr,
                              ampliconNames,
                              minimumMappingQuality = 20,
                              removeDup = FALSE) {

  if (missing(sampleNames)) {
    sampleNames = bamFilenames
  }

  if (missing(ampliconNames)) {
    ampliconNames = elementMetadata(gr)$ampliconNames
  }

  bamIndexFilepaths <- BamIndexFilepaths(bamFilenames)

  # Because of package check complaints
  i <- NULL
  curbam = foreach(i = seq_along(bamFilenames), .combine = cbind) %do% {

    if (!file.exists(bamIndexFilepaths[i])) {
      message(paste("Bai file not found. Creating bai file for", bamFilenames[i], "..."))
      IndexMultipleBams(bamFilenames)
    }

    message(paste("Reading counts for bam file: ", bamFilenames[i]))
    countBamInGRanges(bamFilenames[i],
                      gr,
                      remove.dup = removeDup,
                      min.mapq = minimumMappingQuality,
                      get.width = TRUE)
  }

  listOfBamsHasOnlyOneElement <- length(bamFilenames) == 1
  if (listOfBamsHasOnlyOneElement) {
    curbam <- matrix(curbam)
  }

  colnames(curbam) = sampleNames
  rownames(curbam) = ampliconNames
  return(curbam)
}

#get the combined counts
CombinedNormalizedCounts <- function (sampleCounts,
                                      referenceCounts,
                                      method = "tmm",
                                      ampliconNames = NULL) {
  sampleCounts <- as.matrix(sampleCounts)
  referenceCounts <- as.matrix(referenceCounts)
  allCounts <- cbind(sampleCounts, referenceCounts)
  classes <- rep(c("samples", "reference"), c(ncol(sampleCounts),
                                              ncol(referenceCounts)))
  bamDataRangesNorm <- NormalizeCounts(allCounts, method = method)
  #  bamDataRangesNorm <- NormalizeCounts(allCounts)
  if (!is.null(ampliconNames)) {
    rownames(bamDataRangesNorm) = ampliconNames
  }
  normalizedSamples <- bamDataRangesNorm[, which(classes == "samples"), drop = FALSE]
  normalizedReference <- bamDataRangesNorm[, which(classes == "reference"), drop = FALSE]
  return(list(samples = normalizedSamples, reference = normalizedReference))
}

###############################################################################
# helper functions
###############################################################################

VerifiyIfOutputDirectoryExistsOrIsNotEmpty <- function(outputDir) {
  outputDirectoryExists <- file.exists(outputDir)
  outputDirectoryIsNotEmpty <- length(dir(outputDir, all.files=TRUE)) != 0
  if (outputDirectoryExists || outputDirectoryIsNotEmpty) {
    stop("Output directory already exists or is not empty. Please delete or change output directory")
  }
}

tss <- function(allCounts) {
  return(apply(allCounts, 2, function(x) {x / sum(x)}))
}

NormalizeCounts <- function(allCounts, method = "tmm") {
  # TODO check if normalization method is supported
  allowedNormalizationMethods <- c("tmm", # trimmed mean of M-values
                                   "tss") # total sum scaling
  tinyValueToAvoidZeroReadCounts <- 1e-05
  bamDataRangesNorm <- NULL
  if (method == "tmm") {
    bamDataRangesNorm <- tmm(allCounts, refColumn = NULL, k = tinyValueToAvoidZeroReadCounts)
  } else if (method == "tss") {
    #    bamDataRangesNorm <- apply(allCounts, 2, function(x) {x / sum(x)})
    bamDataRangesNorm <- tss(allCounts)
  } else {
    message(paste("method", method, "not included in the allowed Normalization methods", allowedNormalizationMethods))
  }
  return(bamDataRangesNorm)
}

#get the positions for each gene as a list
IndexGenesPositions = function(genes) {
  positions = seq_along(genes)
  genesPos = split(positions, genes)
  return(genesPos)
}

# #how many numbers in a vector are above a given cutoff
# PercentAboveValue <- function(vector, value) {
#     sum(vector > value)/length(vector)
# }
#
# #how many numbers in a vector are below a given cutoff
# PercentBelowValue <- function(vector, value) {
#     sum(vector < value)/length(vector)
# }

BamIndexFilepaths <- function(bamFilepaths, indexType = ".bam.bai", ignoreCase = FALSE) {
  return (gsub(".bam$",
               bamFilepaths,
               replacement = indexType,
               ignore.case = ignoreCase))
}

# check this params with Thomas
#index the bam files if there is no index yet
IndexMultipleBams <- function(bams,
                              index_type = ".bam.bai") {
  #check if the index already exists and need to be indexed
  #     potentialBaiFilenames <- gsub(".bam$",
  #                                 bams,
  #                                 replacement = index_type,
  #                                 ignore.case = TRUE)

  potentialBaiFilenames <- BamIndexFilepaths(bams, indexType = index_type)

  print(paste("potentialBaiFilenames", potentialBaiFilenames))

  bamsToBeIndexed <- bams[!sapply(potentialBaiFilenames, file.exists)]
  if(length(bamsToBeIndexed) > 0) {
    #index the bams
    #if(multicore) {
    # check with Thomas if we can replace by this..
    #  bplapply(bamsToBeIndexed, indexBam)
    #mclapply(bamsToBeIndexed, indexBam, mc.cores = ncores)
    #} else {
    lapply(bamsToBeIndexed, indexBam)
    #}
  }
}

WriteListToXLSX <- function(listOfDataFrames, multipleFiles = FALSE, outputFolder = file.path(getwd(), "xlsx"), filepath = "list.xlsx") {
  if (multipleFiles) {
    dir.create(outputFolder, recursive = TRUE, showWarnings = FALSE)
    for(name in names(listOfDataFrames)){
      message(paste0("Saving file to '", file.path(outputFolder, paste0(name, ".xlsx"))))
      write.xlsx(listOfDataFrames[[name]], file = file.path(outputFolder, paste0(name, ".xlsx")), rowNames = TRUE)
    }
  } else {
    message(paste0("Saving file to '", filepath, "'"))
    write.xlsx(listOfDataFrames, file = filepath, rowNames = TRUE)
  }
}

ReadXLSXToList <- function(filepath, rowNames = TRUE, colNames = TRUE) {
  listOfDataFrames <- list()
  for(name in getSheetNames(filepath)) {
    listOfDataFrames[[name]] <- read.xlsx(filepath,
                                          rowNames = rowNames,
                                          colNames = colNames,
                                          sheet = name)
  }
  return(listOfDataFrames)
}

###############################################################################
# functions for bootstrapping
###############################################################################

adjustedLength <- function(a) {
  round(sqrt(length(a)))
}

ratiosMean <- function(ratios) {
  return(exp(mean(log(ratios))))
}

#a function that randomly samples positions from a vector using subsampling
SubsamplingPositions <- function(pos, mtry = adjustedLength(pos)) {
  return(sample(pos, mtry, replace = FALSE))
}

#calculates the mean from a vector given specific positions
PositionMean <- function(position, vector) {
  return(ratiosMean(vector[position]))
}

#calculate the mean for each gene
GenePositionMean <- function(genesPos, vector) {
  means = sapply(genesPos, PositionMean, vector = vector)
  names(means) = names(genesPos)
  return(means)
}

#calculates the mean from a vector given specific positions for a matrix
GeneMeanRatioMatrix <- function(genesPos, ratioMatrix) {
  # Package check complaints
  i <- NULL
  ratioList =  foreach(i = 1:ncol(ratioMatrix), .combine = cbind) %do% {
    GenePositionMean(genesPos,ratioMatrix[, i])
  }
  ratioList <- as.matrix(ratioList)
  colnames(ratioList) <- colnames(ratioMatrix)
  return(ratioList)
}

# ReferenceWeights <- function(refmat, varianceFunction = sd) {
#     variance = apply(refmat, 2, varianceFunction)
#     weights = 1/variance
#     return(weights)
# }

# split into a function generating the bootstrap distribution
# and a function averaging over amplicons.
# define the function to generate bootstrap distributions and
# return a list with the genewise bootstrapping
BootList <- function(geneNames, sampleMatrix, refmat, replicates) {
  enforceDeterministicResult()
  # get the genes positions in the matrix as a list from a gene name vector
  genesPos <- IndexGenesPositions(geneNames)

#  if (class(sampleMatrix) != "matrix") {
  if (!("matrix" %in% class(sampleMatrix))) {
        stop(paste("Parameter 'sampleMatrix' has to be of class matrix and is of class", class(sampleMatrix)))
  }

  if (is.null(colnames(sampleMatrix)) | length(colnames(sampleMatrix))!=ncol(sampleMatrix)) {
    stop("All columns of 'sampleMatrix' have to be named")
  }

  # a vector and matrix are not the same and for a vector iterating over
  # the column makes no sense so we have
  # to check if a matrix or a vector was passed. ncol only works for matrix
  # not for vector
#  if (class(sampleMatrix) == "matrix") {
  if ("matrix" %in% class(sampleMatrix)) {
    iterator <- 1:ncol(sampleMatrix)
    # } else {
    #     iterator <- 1
  }

  i <- NULL
  j <- NULL
  bootListSamples <- foreach(i = iterator) %:%
    foreach(j = rep(1, replicates), .combine = rbind) %do% {
      # a vector and matrix are not the same and for a vector
      # iterating over the column makes no sense so we have to check
      # if a mtrix or a vector was passed.
#      if (class(sampleMatrix) == "matrix") {
      if ("matrix" %in% class(sampleMatrix)) {
        testSample <- sampleMatrix[, i]
        # } else {
        #     testSample <- sampleMatrix
      }
      # for each gene subsample the amplicon positions independently
      # sample the samples using bootstrapping
      sampleBootPos <- sample(1:ncol(refmat),
                              ncol(refmat),
                              replace = TRUE)
      geneBootPos <- c(lapply(genesPos,
                              SubsamplingPositions),
                       recursive = TRUE)
      # given the obtained sampling using the bootstraps calculated
      refMatPos <- refmat[geneBootPos, sampleBootPos, drop=FALSE]

      bootRatio <- testSample[geneBootPos]/rowMeans(refMatPos)

      # after the bootstrapping the gene positions in the vector
      # changes so recalculate them
      splitClass <- rep(names(genesPos), sapply(genesPos,
                                                adjustedLength))
      newGenesPos <- split(seq_along(splitClass), splitClass)
      sapply(newGenesPos, PositionMean, vector = bootRatio)
    }

  names(bootListSamples) <- basename(colnames(sampleMatrix))

  # In case a single step of bootstrap is done (for testing purposes),
  # the elements of each sample will be converted to a matrix
  for (i in seq_along(bootListSamples)) {
    bootListSamples[[i]] <- matrix(bootListSamples[[i]], ncol = length(names(genesPos)))
    colnames(bootListSamples[[i]]) <- names(genesPos)
  }
  return(bootListSamples)
}



ConsensusCheck <- function(genomicRangesFromBed, sampleMatrix, referenceMatrix) {
  #  sampleMatrix <- samplesNormalizedReadCounts
  #  referenceMatrix <- referenceNormalizedReadCounts

  #  View(sampleMatrix[, c(1,2)])
  #  View(apply(referenceMatrix, 1, median))

  #  ConsensusCheck(genomicRangesFromBed, results_1PGM_sequencingLibraries1@sampleReadCounts, results_1PGM_sequencingLibraries1@referenceReadColunts),

  #  sampleMatrix <- results_1PGM_sequencingLibraries1@sampleReadCounts
  #  referenceMatrix <- results_1PGM_sequencingLibraries1@referenceReadColunts


  # This should not reach here.. but just in case..
  colnames(sampleMatrix) <- basename(colnames(sampleMatrix))
  colnames(referenceMatrix) <- basename(colnames(referenceMatrix))

  geneNames <- genomicRangesFromBed$geneNames

  allSampleRatios <- data.frame()
  for (smpl in colnames(sampleMatrix)) {
    #    smpl <- colnames(sampleMatrix)[which(grepl("28208", colnames(sampleMatrix)))]
    ampliconsRatio <- sampleMatrix[, smpl]/apply(referenceMatrix, 1, median)
    allSampleRatios <- rbind(allSampleRatios, ampliconsRatio)
  }

  allSampleRatios <- t(allSampleRatios)
  rownames(allSampleRatios) <- rownames(sampleMatrix)
  colnames(allSampleRatios) <- colnames(sampleMatrix)

  genesPositions <- IndexGenesPositions(geneNames)

  geneConsensus <- data.frame(matrix(nrow=length(genesPositions), ncol=ncol(sampleMatrix)))
  geneConsensus <- data.frame()
  for (smpl in colnames(allSampleRatios)) {
    for (gene in names(genesPositions)) {
      geneAmpliconsRatios <- allSampleRatios[,smpl][genesPositions[[gene]]]
      ampliconConsensus <- all(geneAmpliconsRatios>1) | all(geneAmpliconsRatios<1)
      geneConsensus[gene, smpl] <- ampliconConsensus
    }
  }
  # TODO case the reference has amplicons with 0 read counts, it returns NA instead of boolean
  geneConsensus[is.na(geneConsensus)] <- FALSE
  return(geneConsensus)
}

HaveMininumNumberOfAmplicons <- function(genomicRangesFromBed, numberOfAmplicons) {
  geneNames <- genomicRangesFromBed$geneNames
  genesPositions <- IndexGenesPositions(geneNames)
  sapply(genesPositions, function(x) {length(x)>=numberOfAmplicons})
}

#calculate significance
CheckSignificance <- function(bootList, significanceLevel = 0.05) {
  # Bonferroni correction
  significanceLevel <- significanceLevel/ncol(bootList[[1]])
  margin <- significanceLevel/2

  # package check complains
  i <- NULL
  j <- NULL
  sigTables =  foreach(i = seq_along(bootList)) %:%
    foreach(j = 1:ncol(bootList[[i]]), .combine = rbind) %do% {
      bootRatioDistribution <- bootList[[i]][,j]
      #      meanNoise <- exp(mean(log(bootRatioDistribution)))
      meanNoise <- ratiosMean(bootRatioDistribution)
      lowerBound <- quantile(bootRatioDistribution, margin, type = 1)
      upperBound <- quantile(bootRatioDistribution, 1 - margin, type = 1)
      bounds <- c(lowerBound, meanNoise, upperBound)
      roundedBounds <- round(bounds, digits = 2)
      maybeAmplification <- roundedBounds[1] > 1
      maybeDeletion <- roundedBounds[3] < 1
      roundedSignificant <- maybeAmplification | maybeDeletion
      copyNumberPutativeStatus <- "Normal"
      if(maybeAmplification) {
        copyNumberPutativeStatus <- "Amplification"
      } else {
        if (maybeDeletion) {
          copyNumberPutativeStatus <- "Deletion"
        }
      }
      sigTable <- c(roundedBounds, roundedSignificant, copyNumberPutativeStatus)
      names(sigTable)  <- c("lowerBound", "mean", "upperBound", "isSig", "putativeStatus")
      sigTable
    }
  names(sigTables) <- names(bootList)

  # name the genes in the list. All samples should share the same gene names
  geneNames <- colnames(bootList[[1]])
  for (i in seq_along(sigTables)) {
    rownames(sigTables[[i]]) <- geneNames
  }
  return(sigTables)
}

#SignificantGenes <- function(sigList, genesPositionsIndex) {
#  # package check complains
#  i <- NULL
#  sigGenes = foreach(i = seq_along(sigList)) %do% {
#    names(genesPositionsIndex)[sigList[[i]][, "isSig"] == TRUE]
#  }
#  return(sigGenes)
#}

# #############################################################################
# # functions for background noise estimation
# #############################################################################
NonSignificantGeneIndex <- function(sigList, genesPositionsIndex) {
  # package check complains
  i <- NULL
  genePosNonSig = foreach(i = seq_along(sigList)) %do% {
    sigGenes = names(genesPositionsIndex)[sigList[[i]][, "isSig"] == TRUE]
    selectedIndex <- RemSigGenes(genesPositionsIndex, sigGenes)
    if (length(selectedIndex) == 0) {
      genesPositionsIndex
    } else {
      selectedIndex
    }
  }
  names(genePosNonSig) <- names(sigList)
  return(genePosNonSig)
}

AmplProbMultipeSamples <- function(genePosNonSig) {
  # package check complains
  i <- NULL
  amplWeights = foreach(i = seq_along(genePosNonSig)) %do% {
    AmplProb(genePosNonSig[[i]])
  }
  names(amplWeights) <- names(genePosNonSig)
  return(amplWeights)
}

enforceDeterministicResult <- function() {
  set.seed(1)
}

Background <- function(geneNames,
                       samplesNormalizedReadCounts,
                       referenceNormalizedReadCounts,
                       bootList,
                       replicates = 1000,
                       significanceLevel = 0.05,
                       robust = FALSE) {
  enforceDeterministicResult()
  #which genes showed significant changes
  sigList <- CheckSignificance(bootList)
  # gene index
  genesPositionsIndex <- IndexGenesPositions(geneNames)
  # calculate the reference mean
  #refMean <- apply(referenceNormalizedReadCounts, 1, mean)
  refMean <- rowMeans(referenceNormalizedReadCounts)
  # calculate the ratio matrix for each sample
  ratioMatrix <- RatioMatrix(samplesNormalizedReadCounts, refMean)
  # remove the significant genes from the noise estimation
  genesPosNonSig <- NonSignificantGeneIndex(sigList, genesPositionsIndex)
  # calculate the weight for each amplicon of non significant changed genes
  amplWeights <- AmplProbMultipeSamples(genesPosNonSig)
  uniqueAmpliconNumbers <- NumberOfUniqueAmplicons(genesPositionsIndex)

  sampleIndex <- NULL
  backgroundObject <- foreach(sampleIndex = seq_along(genesPosNonSig)) %do% {
    message(paste0("Calculating Background for ", colnames(samplesNormalizedReadCounts)[sampleIndex]))
    nonSigAmpliconRatios <- ratioMatrix[unlist(genesPosNonSig[[sampleIndex]]), sampleIndex]
    IterateAmplNum(uniqueAmpliconNumbers,
                   nonSigAmpliconRatios,
                   replicates = replicates,
                   probs = amplWeights[[sampleIndex]],
                   significanceLevel = significanceLevel,
                   robust = robust)
  }
  return(backgroundObject)
}

roundDf <- function(x, digits) {
  # round all numeric variables
  # x: data frame
  # digits: number of digits to round
  #  numeric_columns <- sapply(x, mode) == 'numeric'
  numeric_columns <- sapply(x, is.numeric)
  x[numeric_columns] <-  round(x[numeric_columns], digits)
  x
}

# Helper function to add column to reportTables
ReliableStatus <- function(putative, numberOfThresholdsPassed) {
  status <- as.vector(putative)
  #  print(status)
  if (numberOfThresholdsPassed < 2) {
    status <- "Normal"
  }
  return(status)
}

ReportTables <- function(geneNames,
                         samplesNormalizedReadCounts,
                         referenceNormalizedReadCounts,
                         bootList,
                         backgroundNoise) {
  # get the background noise in a format that can be used
  # for a report table
  backgroundReport <- BackgroundReport(backgroundNoise, geneNames)
  # gene index
  genesPositionsIndex <- IndexGenesPositions(geneNames)
  # calculate the reference mean
  #refMean <- apply(referenceNormalizedReadCounts, 1, mean)
  refMean <- rowMeans(referenceNormalizedReadCounts)
  # calculate the ratio matrix for each sample
  ratioMatrix <- RatioMatrix(samplesNormalizedReadCounts, refMean)
  # calculate the genewise ratio matrix from the ratio_mat
  ratioMatGene <- GeneMeanRatioMatrix(genesPositionsIndex, ratioMatrix)
  sigList <- CheckSignificance(bootList)
  # because package generaton complains..
  i <- NULL
  reportTables <- foreach(i = seq_along(backgroundReport)) %do% {
    #isSig <- (sigList[[i]][, "upperBound"] < 1) | (sigList[[i]][, "lowerBound"] > 1)
    putativeStatus <- sigList[[i]][, "putativeStatus"]
    isSig <- putativeStatus != "Normal"

    backgroundUp <- backgroundReport[[i]][, "UpperNoise"]
    backgroundDown <- backgroundReport[[i]][, "LowerNoise"]
    rMatGene <- ratioMatGene[, i]
    # TODO is this correct?! why compare rMatGene > 1 ?
    aboveNoise <- (rMatGene > 1 & sigList[[i]][, "lowerBound"] > backgroundUp) |
      (rMatGene < 1 & sigList[[i]][, "upperBound"] < backgroundDown)

    dfTemp <- data.frame(meanRatio = ratioMatGene[, i],
                         lowerBoundBootstrapRatio = sigList[[i]][, "lowerBound"],
                         meanBootstrapRatio = sigList[[i]][, "mean"],
                         upperBoundBootstrapRatio = sigList[[i]][, "upperBound"],
                         #                          lowerNoise = backgroundReport[[i]][, 1],
                         #                          meanNoise = backgroundReport[[i]][, 2],
                         #                          upperNoise = backgroundReport[[i]][, 3],
                         lowerNoise = backgroundReport[[i]][, "LowerNoise"],
                         meanNoise = backgroundReport[[i]][, "MeanNoise"],
                         upperNoise = backgroundReport[[i]][, "UpperNoise"],
                         significant = isSig,
                         aboveNoise = aboveNoise,
                         amplNum = as.vector(table(geneNames)),
                         putativeStatus = putativeStatus,
                         passed = isSig + aboveNoise)

    significativeNumbers <- 2
    dfTemp <- roundDf(dfTemp, significativeNumbers)
    names(dfTemp) <- c("MeanRatio",
                       "LowerBoundBoot",
                       "MeanBoot",
                       "UpperBoundBoot",
                       "LowerNoise",
                       "MeanNoise",
                       "UpperNoise",
                       "Signif.",
                       "AboveNoise",
                       "Amplicons",
                       "PutativeStatus",
                       "Passed")
    dfTemp
  }

  # TODO Check if this column added is correct
  for(i in 1:length(reportTables)) {
    tmpReportTable <- reportTables[[i]]
    #  tmpReportTable <- reportTables[[1]]
    status <- mapply(ReliableStatus, tmpReportTable$PutativeStatus, tmpReportTable$Passed)
    #  print(i)
    reportTables[[i]] <- cbind(reportTables[[i]], "ReliableStatus" = status)
  }

  #    names(reportTables) <- colnames(samplesNormalizedReadCounts)
  names(reportTables) <- names(bootList)
  return(reportTables)
}

BackgroundReport <- function(background, geneNames) {
  # because package generation complains..
  i <- NULL
  j <- NULL
  backgroundReport = foreach(i = seq_along(background)) %:%
    foreach(j = table(geneNames), .combine = rbind) %do% {
      background[[i]][[as.character(j)]]
    }
}

#remove genes that were considered significant by the algorithm
RemSigGenes <- function(genesPos, sigGenes) {
  #if some of the genes were reported with a significantly different read count
  #should not be used for
  #background testing
  #which genes are in the genes pos object
  geneNames= names(genesPos)
  #define a function for not in
  `%ni%` = Negate(`%in%`)
  #which of these genes did not show a significant read count  change
  nonSigGenes = which(geneNames %ni% sigGenes)
  #keep only the genes that did not show a significant read coount change
  nonSigGenesPos = genesPos[nonSigGenes]
  return(nonSigGenesPos)
}

AmplProb <- function(genesPos) {
  #how many amplicons where used for each ofhte genes
  geneCounts = countAmplicons(genesPos)
  #adjust the probablity for depending on the number of amplicons for each gene
  genePerc =  1/geneCounts
  #this information has to be available for each stable position
  ampliconProb = rep(genePerc, geneCounts)
  return(ampliconProb)
}

RatioMatrix <- function(sampleMat, refMean) {
  apply(sampleMat, 2, `/`, refMean)
}

SampleRatio <- function(ratios, numAmpl, amplWeights = NULL) {
  replace <- FALSE
  if(numAmpl >= length(ratios)) {
    replace <- TRUE
  }
  randomPos = sample(seq_along(ratios),
                     numAmpl,
                     prob = amplWeights,
                     replace = replace)
  randomlySelectedRatios = ratios[randomPos]
  #  randomMean = exp(mean(log(randomlySelectedRatios)))
  randomMean <- ratiosMean(randomlySelectedRatios)
  return(randomMean)
}

DescriptiveStatistics <- function(distribution, robust = FALSE) {
  #	print(paste("Robust is ", robust))
  if (robust) {
    centralTendency <- median(distribution)
    variability <- mad(distribution, constant = 1) # median absolute deviation
  } else {
    centralTendency <- mean(distribution)
    variability <- sd(distribution)
  }
  return(list(centralTendency = centralTendency,
              variability = variability))
}

SampleNoiseGenes <- function(numAmpl = 2,
                             ratios,
                             replicates = 100,
                             probs = NULL,
                             significanceLevel = 0.05,
                             robust = FALSE) {
  margin <- significanceLevel / 2
  # now repeat the sampling for the selected number of
  # amplicons replicates
  sampleNoiseDistribution = replicate(replicates,
                                      SampleRatio(ratios,
                                                  numAmpl,
                                                  amplWeights = probs))
  #get the upper, mean and lower values of the noise distribution
  logSampleNoiseDistribution <- log(sampleNoiseDistribution)

  ds <- DescriptiveStatistics(logSampleNoiseDistribution, robust)
  logMeanNoise <- ds$centralTendency
  logSdNoise <- ds$variability

  # if (robust) {
  #     logMeanNoise <- median(logSampleNoiseDistribution)
  #     logSdNoise <- mad(logSampleNoiseDistribution, constant = 1)
  # } else {
  #     logMeanNoise <- mean(logSampleNoiseDistribution)
  #     logSdNoise <- sd(logSampleNoiseDistribution)
  # }

  lowerBound <- exp(logMeanNoise + qnorm(margin) * logSdNoise)
  meanNoise <- exp(logMeanNoise)
  upperBound <- exp(logMeanNoise + qnorm(1 - margin) * logSdNoise)

  sampledNoise = c(lowerBound, meanNoise, upperBound)
  names(sampledNoise) = paste0(c("Lower", "Mean", "Upper"), "Noise")
  return(sampledNoise)
}

countAmplicons <- function(x) {
  return(sapply(x, NROW))
}

# Returns a list with the unique number of amplicons for all genes
NumberOfUniqueAmplicons <- function(genesPos) {
  # calcUniqueAmpliconNumbers <-function(genesPos) {
  ampliconNumbers = countAmplicons(genesPos)
  uniqueAmpliconNumbers = sort(unique(ampliconNumbers))
  return(uniqueAmpliconNumbers)
}

#now calculate the background for each of the unique amplicon numbers
# ratios: the
IterateAmplNum <- function(uniqueAmpliconNumbers,
                           ratios,
                           replicates = 100,
                           probs = NULL,
                           significanceLevel = 0.05,
                           robust = FALSE) {
  # Needed because of package check complaints..
  i <- NULL
  noiseResults = foreach(i = seq_along(uniqueAmpliconNumbers)) %do% {
    sampledNoise = SampleNoiseGenes(uniqueAmpliconNumbers[i],
                                    ratios = ratios,
                                    replicates = replicates,
                                    probs = probs,
                                    significanceLevel = significanceLevel,
                                    robust = robust)
    sampledNoise
  }
  names(noiseResults) = as.character(uniqueAmpliconNumbers)
  return(noiseResults)
}

SelectReferenceSetByPercentil <- function(allSamplesReadCounts,
                                          normalizationMethod = "tmm",
                                          lowerBoundPercentage = 1,
                                          upperBoundPercentage = 99) {

  allSamplesReadCountsNormalized <-NormalizeCounts(allSamplesReadCounts,
                                                   method = normalizationMethod)

  selectedSamplesFilepath <- NULL
  allSamples <- round(allSamplesReadCountsNormalized)
  samplesColnames <- colnames(allSamples)

  selectedRange <- list()
  for (i in 1:nrow(allSamples)) {
    #    selectedRange[[i]] <- quantile(allSamples[i,])
    selectedRange[[i]] <- quantile(allSamples[i,], c(lowerBoundPercentage, upperBoundPercentage) / 100)
  }

  selectedSamplesIndex <- NULL
  excludedSamplesIndex <- NULL
  for (j in 1:ncol(allSamples)) {
    addSample <- TRUE
    for(i in 1:nrow(allSamples)) {
      if ((allSamples[i, j] < selectedRange[[i]][1]) | (allSamples[i, j] > selectedRange[[i]][2])) {
        addSample <- FALSE
        #        print(paste("Add sample?", addSample, i, selectedRange[[i]][1], allSamples[i, j], selectedRange[[i]][2]))
        break
      }
    }
    if (addSample) {
      selectedSamplesIndex <- cbind(selectedSamplesIndex, j)
    } else {
      excludedSamplesIndex <- cbind(excludedSamplesIndex, j)
    }
  }

  selectedSamples <- allSamples[, selectedSamplesIndex]
  excludedSamples <- allSamples[, -excludedSamplesIndex]
  selectedSamplesFilepath <- colnames(selectedSamples)
  return(selectedSamplesFilepath)
}

SelectReferenceSetByInterquartileRange <- function(allSamplesReadCounts, normalizationMethod = "tmm", iqrFactor = 1) {
  allSamplesReadCounts <- NormalizeCounts(allSamplesReadCounts, method = normalizationMethod)
  fivenumReads <- t(apply(allSamplesReadCounts, 1, fivenum))
  colnames(fivenumReads) <- c("SampleMinimum", "LowerQuartile", "Median", "UpperQuartile", "SampleMaximum")
  iqr <- fivenumReads[, "UpperQuartile"] - fivenumReads[, "LowerQuartile"]
  dispersion <- iqrFactor * iqr
  sampleIndexes <- NULL
  for (i in 1:ncol(allSamplesReadCounts)) {
    if (all((allSamplesReadCounts[, i] >= (fivenumReads[,"LowerQuartile"] - dispersion ))  & ((dispersion + fivenumReads[,"UpperQuartile"])  >= allSamplesReadCounts[, i]))) {
      sampleIndexes <- c(sampleIndexes, i)
    }
  }
  return(colnames(allSamplesReadCounts)[sampleIndexes])
}

SelectReferenceSetByKmeans <- function(allSamplesReadCounts, normalizationMethod = "tmm", referenceNumberOfElements) {
  # referenceSamples <- function(referenceNumberOfElements, allSamples) {

  allSamples <- NormalizeCounts(allSamplesReadCounts, method = normalizationMethod)

  #  allSamples <- allSamplesReadCounts
  numberOfSamples <- ncol(allSamples)
  numberOfOutliers <- numberOfSamples - referenceNumberOfElements
  if (numberOfOutliers < 0) {
    message(paste("Number of samples", numberOfSamples, "is smaller than the size of the reference set requested", referenceNumberOfElements, ", so all samples will be selected"))
    return(colnames(allSamples))
  } else if (numberOfOutliers == 0) {
    message(paste("Number of samples", numberOfSamples, "is equal to the size of the reference set requested", referenceNumberOfElements))
    return(colnames(allSamples))
  } else {
    allSamples <- t(allSamples)
    kmeans.result <- kmeans(allSamples, centers=1)
    # "centers" is a data frame with one center
    centers <- kmeans.result$centers[kmeans.result$cluster, ]
    distances <- sqrt(rowSums((allSamples - centers)^2))
    outliers <- order(distances, decreasing=TRUE)[1:numberOfOutliers]
    # these rows are top outliers
    message(cat("outliers:", colnames(allSamplesReadCounts)[outliers]))
    allSamples <- t(allSamples)
    #    return(colnames(allSamples[, -outliers, drop = FALSE]))
    return(colnames(allSamples)[-outliers])
  }
}

SelectReferenceSet <- function(samplesFilepaths,
                               genomicRanges,
                               removePcrDuplicates = FALSE,
                               normalizationMethod = "tmm",
                               referenceMaximumNumberOfElements = 10,
                               referenceSelectionMethod = "percentil",
                               lowerBoundPercentage = 1,
                               upperBoundPercentage = 99) {

  #  samplesFilepaths <- allSamplesFilteredFilepaths[1:3]
  #  genomicRanges <- genomicRangesFromBed

  allSamplesReadCounts <- ReadCountsFromBam(bamFilenames = samplesFilepaths,
                                            sampleNames = samplesFilepaths,
                                            gr = genomicRanges,
                                            #                                            ampliconNames = ampliconNames,
                                            removeDup = removePcrDuplicates)

  selectedSamplesFilepath <- SelectReferenceSetFromReadCounts(allSamplesReadCounts,
                                                              #                                                              genomicRanges,
                                                              #                                                              removePcrDuplicates = removePcrDuplicates,
                                                              normalizationMethod = normalizationMethod,
                                                              referenceMaximumNumberOfElements = referenceMaximumNumberOfElements,
                                                              referenceSelectionMethod = referenceSelectionMethod,
                                                              lowerBoundPercentage = lowerBoundPercentage,
                                                              upperBoundPercentage = lowerBoundPercentage)
  return(selectedSamplesFilepath)
}

SelectReferenceSetFromReadCounts <- function(allSamplesReadCounts,
                                             #                               genomicRanges,
                                             #                               removePcrDuplicates = FALSE,
                                             normalizationMethod = "tmm",
                                             referenceMaximumNumberOfElements = 30,
                                             #                               referenceSelectionMethod = "percentil",
                                             referenceSelectionMethod = "kmeans",
                                             lowerBoundPercentage = 1,
                                             upperBoundPercentage = 99) {


  selectedSamplesFilepath <- NULL
  if (referenceSelectionMethod == "percentil") {
    message("Selecting reference by percentil")
    selectedSamplesFilepath <- SelectReferenceSetByPercentil(allSamplesReadCounts,
                                                             normalizationMethod = normalizationMethod,
                                                             #                                                              referenceSelectionMethod = referenceSelectionMethod,
                                                             lowerBoundPercentage = lowerBoundPercentage,
                                                             upperBoundPercentage = upperBoundPercentage)
  } else if (referenceSelectionMethod == "kmeans") {
    message("Selecting reference by kmeans")
    selectedSamplesFilepath <- SelectReferenceSetByKmeans(allSamplesReadCounts,
                                                          normalizationMethod = normalizationMethod,
                                                          referenceMaximumNumberOfElements)
  } else {
    message(paste("method", referenceSelectionMethod, "not supported"))
  }
  return(selectedSamplesFilepath)
}

CollectColumnFromAllReportTables <- function(reportTables, columnName) {
  compilation <- NULL
  mySampleNames <- names(reportTables)
  for(name in mySampleNames) {
    tmpReportTable <- reportTables[[name]]
    compilation <- cbind(compilation, as.vector(tmpReportTable[, columnName]))
  }
  colnames(compilation) <- mySampleNames
  rownames(compilation) <- rownames(reportTables[[1]])
  return(compilation)
}

revalueDF <- function(myDataFrame, mappings) {
  myDataFrameColnames <- colnames(myDataFrame)
  for(columnName in myDataFrameColnames) {
    myDataFrame[, columnName] <- revalue(myDataFrame[, columnName], mappings)
  }
  mode(myDataFrame) <- "numeric"
  return(myDataFrame)
}

## TODO check this function
# ReliableAberrationStatusHeatMap <- function(meanBootResults,
#                                             realiableResults,
#                                             filepath = NULL) {
#   meanBootResults <- melanomaMeanBootResults
#   reliableResults <- melanomaCNVPanelizerResults
#   test_that("testing correct dimensions", {
#     expect_equal(colnames(meanBootResults), colnames(reliableResults))
#     expect_equal(rownames(meanBootResults), rownames(reliableResults))
#   })
#
#   kNormalDefaultValue <- 1.0
#   kAmplificationMaximalValue <- 4.0
#   kAmplificationMinimumValue <- 2.0
#
#   for(i in rownames(reliableResults)) {
#     for(j in colnames(reliableResults)) {
#       if (reliableResults[i,j] == "Normal") {
#         #        print(paste("Normal", i, reliableResults[i,j], "setting ", meanBootResults[i,j],"as", kNormalDefaultValue))
#         meanBootResults[i,j] <- kNormalDefaultValue
#       }
#
#       if ((reliableResults[i,j] == "Amplification") & (meanBootResults[i,j] > kAmplificationMaximalValue)) {
#         print(paste("Too high ", i, reliableResults[i,j], "setting ", meanBootResults[i,j],"as ", kAmplificationMaximalValue))
#         meanBootResults[i,j] <- kAmplificationMaximalValue
#       }
#
#       if ((reliableResults[i,j] == "Amplification") & (meanBootResults[i,j] < kAmplificationMinimumValue)) {
#         print(paste("Amplification too low", i, reliableResults[i,j], "setting ", meanBootResults[i,j],"as", kNormalDefaultValue))
#         meanBootResults[i,j] <- kNormalDefaultValue
#       }
#     }
#   }
#
#   df.team_data <- melt(meanBootResults)
#   colnames(df.team_data) <- c("Gene", "Sample", "Bootratio")
#   g <- ggplot(data = df.team_data,
#               aes(x = Gene, y = Sample)) +
#     #    geom_tile(aes(fill = scale(Bootratio))) +
#     geom_tile(aes(fill = log(Bootratio))) +
#     #      + opts(theme(axis.text.x = element_text(angle = 180, hjust = 1)))
#     #       scale_fill_gradient2(low="darkblue", high="darkgreen", guide="colorbar") +
#     #      scale_fill_gradient2(low="green4", high="red4", guide="colorbar") +
#     scale_fill_gradient2(low = "darkgreen", mid = "white", high = "darkred") +
#     theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
#     theme(axis.text.y = element_text(size=5))
#   if (!missing(filepath)) {
#     ggsave(filepath)
#   }
# }



# why these colors by default:
StatusHeatmap <- function(dfData,
                          statusColors = c("Deletion" = "blue",
                                           "Normal" = "green",
                                           #                                         " " = "green",
                                           "Amplification" = "red"),
                          header = "Status Heatmap",
                          filepath = "CNVPanelizerHeatMap.png") {

  mappings <- seq(length(statusColors))
  #  mappings <- c(-1, 0, 1)
  names(mappings) <- names(statusColors)
  matrixData <- revalueDF(dfData, mappings)

  valuesUsed <- table(matrixData)
  indexOfColorsUsed <- as.numeric(names(valuesUsed))
  colorsRequired <- unname(statusColors[indexOfColorsUsed])

  # TODO find out why do I need this horrible hack when only one status is available
  if (length(colorsRequired) <= 1) {
    colorsRequired <- c("green", "yellow")
  }

  #  colorsRequired <- as.vector(statusColors[as.numeric(names(table(matrixData)))])

  png(filename = file.path(filepath),    # create PNG for the heat map
      width = 5*300,        # 5 x 300 pixels
      height = 5*300,
      res = 300,            # 300 pixels per inch
      pointsize = 8)        # smaller font size

  lmat <- rbind( c(5,3,4), c(2,1,4) )
  lhei <- c(1.5, 4)
  lhei <- c(1, 10)
  lwid <- c(1.5, 4, 0.75)

  myCexRow <- 0.2

  heatmap.2(matrixData,
            main = header,
            # plot layout
            lmat = lmat,
            lhei = lhei,
            lwid = lwid,
            #          scale = "row",
            #          ColSideColors=heatColors,
            key = FALSE, # whether a color-key should be shown
            #          key = TRUE, # whether a color-key should be shown
            #          key.par=list(match=c(0),
            #             amplification_deletion=c(1, 2, 3, 4),
            #             deletion_amplification=c(5, 6)),
            dendrogram = "none",
            tracecol=NA,
            #          density.info=c("histogram","density","none"),
            #            col = unname(statusColors),
            #            col = c("green"),
            col = colorsRequired,    # if the values are not all revalued the colors have to be passed selectively

            #          col = heatColors,
            #          col = heatColors[matchPoints],
            colsep=1:ncol(matrixData),
            rowsep=1:nrow(matrixData),
            margins = c(15, 5),
            Rowv=FALSE,
            Colv=FALSE,
            #            cexRow=myCexRow
  )

  #   legend("topleft", inset=.02, title= "UOUOU",
  #          "this is legend",
  #          #         fill=topo.colors(3),
  # #         fill=unname(statusColors),
  #          horiz=TRUE, cex=0.8)
  #bottomleft
  #  legend("bottom", inset=.01, title="CNV Status",





  #
  #
  #          legend("left", inset=.01, title="CNV Status",
  # #                legend("left", inset=0, title="CNV Status",
  #          names(statusColors),
  #          #         fill=topo.colors(3),
  #          fill=unname(statusColors),
  #          cex=0.8)
  # #         horiz=TRUE, cex=0.8)
  #




  #   legend("bottomleft", legend=c("Line 1", "Line 2"),
  #          col=c("red", "blue"), lty=1:2, cex=0.8)

  #   legend(1, 95, legend=c("Line 1", "Line 2"),
  #          col=c("red", "blue"), lty=1:2, cex=0.8)

  # in case of mistake this instruction should be called until return null instead of png..
  dev.off()
}

## TODO check
# StatusStability <- function(geneNames, sampleNormalizedReadCounts, tmpReferenceNormalizedReadCounts) {
#   centralTendency <-  matrix(rowMeans(tmpReferenceNormalizedReadCounts), ncol = 1)
#   rownames(centralTendency) <- rownames(tmpReferenceNormalizedReadCounts)
#   referenceAndSampleDifferences <- centralTendency - sampleNormalizedReadCounts
#   rownames(referenceAndSampleDifferences) <- rownames(tmpReferenceNormalizedReadCounts)
#   colnames(referenceAndSampleDifferences) <- "difference"
#   uniqueGeneNames <- unique(geneNames)
#   geneStabilityStatus <- c()
#   for (i in seq(uniqueGeneNames)) {
#     geneAmpliconIndexes <- grep(uniqueGeneNames[i], rownames(myReference))
#     geneAmpliconStatus <- referenceAndSampleDifferences[geneAmpliconIndexes, ]
#     allGeneAmpliconsHaveSameStatus <- length(unique(geneAmpliconStatus >= 0))==1
#     geneStabilityStatus <- c(geneStabilityStatus, allGeneAmpliconsHaveSameStatus)
#   }
#   names(geneStabilityStatus) <- uniqueGeneNames
#   return(geneStabilityStatus)
# }

## TODO CHECK
# StatusStabilityTable <- function(geneNames, tmpSamplesNormalizedReadCounts, tmpReferenceNormalizedReadCounts) {
#   statusStabilityTable <- NULL
#   for (i in seq(ncol(tmpSamplesNormalizedReadCounts))) {
#     statusStabilityTable <- cbind(statusStabilityTable, StatusStability(geneNames, matrix(tmpSamplesNormalizedReadCounts[, i], ncol = 1), tmpReferenceNormalizedReadCounts))
#   }
#   colnames(statusStabilityTable) <- colnames(tmpSamplesNormalizedReadCounts)
#   #   statusStabilityTable <- t(statusStabilityTable)
#   #   colnames(statusStabilityTable) <- colnames(tmpSamplesNormalizedReadCounts)
#   return(statusStabilityTable)
# }

PlotBootstrapDistributions  <- function(bootList,
                                        reportTables,
                                        outputFolder = getwd(),
                                        sampleNames = NULL,
                                        save = FALSE,
                                        #                                        scale = 7) {
                                        scale = 10) {
  selSample <- NULL
  plotList <- foreach(selSample = seq_along(bootList)) %do% {
    test <- as.factor(reportTables[[selSample]][, "Passed"])

    levelLabels <- c("NoChange", "NonReliableChange", "ReliableChange")
    names(levelLabels) <- c(0, 1, 2)
    test <- suppressMessages(revalue(test, levelLabels))
    namedColors <- c("#56B4E9", "#CC79A7" ,"#D55E00")
    names(namedColors) <- levelLabels

    ratios <- NULL
    testsPassed <- NULL
    df <- data.frame(class = as.factor(colnames(bootList[[selSample]])),
                     ratios = (as.vector(t(bootList[[selSample]]))),
                     testsPassed = test)

    #if a genomic ranges object has been
    #if(!is.null(gr)) {
    # newGeneOrder <- bedToGeneOrder(gr)
    # df$class = with(df,factor(class,levels(class)[newGeneOrder]))
    #}

    ylim1 <- boxplot.stats((df$ratios))$stats[c(1, 5)]
    ylim1 <- c(max(ylim1),max(ylim1))
    ylim1 <- log2(ylim1)
    ylim1 <- ylim1 * c(-scale,scale)

    if(is.null(sampleNames)) {
      filename <- names(bootList[selSample])
      filepath <- paste0(outputFolder, "/", filename, "_plot.pdf")
    } else {
      filename <- sampleNames[selSample]
      filepath <- paste0(outputFolder, "/", filename, ".pdf")
    }
    bootPlot <- ggplot(df, aes(x = class, y = log2(ratios), fill = testsPassed)) +
      geom_violin() + ggtitle(filename) +
      theme(plot.title = element_text(lineheight = 0.8, face = "bold"),
            text = element_text(size = 15),
            axis.text.x = element_text(angle = 90)) +
      scale_fill_manual(name = "CNV Reliability",
                        values = namedColors,
                        labels = c("0" = "Foo", "1" = "Bar")) +
      coord_cartesian(ylim = ylim1) +
      theme(axis.text=element_text(size=10),
            axis.title=element_text(size=20,face="bold")) +
      scale_x_discrete("Gene Names") +
      #       geom_hline(yintercept=log2(1.5), color="#009E73") +
      #       geom_hline(yintercept=log2(0.5), color="#009E73") +
      geom_hline(yintercept=0, color="#009E73")

    if(save == TRUE) {
      message(paste0("Saving plot to '", filepath, "'"))
      #            dir.create(outputFolder, recursive = TRUE, showWarnings = TRUE)
      dir.create(outputFolder, recursive = TRUE, showWarnings = FALSE)
      ggsave(filename = filepath,
             plot = bootPlot,
             height = 7.42 * 1,
             width = 8.11 * 2,
             limitsize = FALSE)
    }
    #        dev.off()
    return(bootPlot)
  }
  names(plotList) <- names(reportTables)
  return(plotList)
}

CNVPanelizerFromReadCountsHELPER <- function(sampleReadCounts,
                                             referenceReadCounts,
                                             genomicRangesFromBed,
                                             numberOfBootstrapReplicates = 10000,
                                             normalizationMethod = "tmm",
                                             robust = TRUE,
                                             backgroundSignificanceLevel = 0.05,
                                             outputDir = file.path(getwd(), "CNVPanelizer"),
                                             splitSize = 5) {
  myCNVPanelizerResults <- list()

  splits <- split(1:ncol(sampleReadCounts), ceiling(seq_along(1:ncol(sampleReadCounts))/splitSize))

  #    for (i in 1:ncol(sampleReadCounts)) {

  for (i in 1:length(splits)) {
    myCNVPanelizerResults[[i]] <- CNVPanelizerFromReadCounts(sampleReadCounts[, splits[[i]]],
                                                             referenceReadCounts,
                                                             genomicRangesFromBed,
                                                             numberOfBootstrapReplicates = numberOfBootstrapReplicates,
                                                             normalizationMethod = normalizationMethod,
                                                             robust = robust,
                                                             backgroundSignificanceLevel = backgroundSignificanceLevel,
                                                             outputDir = outputDir)
  }

  n <- length(myCNVPanelizerResults)
  res <- NULL
  for (i in seq(n)) {
    res <- cbind(res, myCNVPanelizerResults[[i]])
  }

  myCNVPanelizerResultsCompiled <- myCNVPanelizerResults[[1]]

  if(length(myCNVPanelizerResults) > 1) {
    for (i in 2:length(myCNVPanelizerResults)) {
      myCNVPanelizerResultsCompiled$sampleReadCounts <- cbind(myCNVPanelizerResultsCompiled$sampleReadCounts, myCNVPanelizerResults[[i]]$sampleReadCounts)
      myCNVPanelizerResultsCompiled$bootList <- append(myCNVPanelizerResultsCompiled$bootList, myCNVPanelizerResults[[i]]$bootList)
      myCNVPanelizerResultsCompiled$backgroundNoise <- append(myCNVPanelizerResultsCompiled$backgroundNoise, myCNVPanelizerResults[[i]]$backgroundNoise)
      myCNVPanelizerResultsCompiled$reportTables <- append(myCNVPanelizerResultsCompiled$reportTables, myCNVPanelizerResults[[i]]$reportTables)
    }
  }
  return(myCNVPanelizerResultsCompiled)
}

######################################################################################
# still on testing..
######################################################################################
# CNVPanelizerFromReadCountsHelperParallel <- function(sampleReadCounts,
#                                                      referenceReadCounts,
#                                                      genomicRangesFromBed,
#                                                      numberOfBootstrapReplicates = 10000,
#                                                      normalizationMethod = "tmm",
#                                                      robust = TRUE,
#                                                      backgroundSignificanceLevel = 0.05,
#                                                      outputDir = file.path(getwd(), "CNVPanelizer"),
#                                                      splitSize = 5) {
#
#   myCNVPanelizerResults <- list()
#
#   splits <- split(1:ncol(sampleReadCounts), ceiling(seq_along(1:ncol(sampleReadCounts))/splitSize))
#
#   # Use the detectCores() function to find the number of cores in system
#   no_cores <- detectCores()
#   # Setup cluster
#   clust <- makeCluster(no_cores - 1)
#   clusterExport(clust,
#                 varlist = c(
#                   "sampleReadCounts",
#                   "referenceReadCounts",
#                   "genomicRangesFromBed",
#                   "numberOfBootstrapReplicates",
#                   "normalizationMethod",
#                   "robust",
#                   "backgroundSignificanceLevel",
#                   "outputDir",
#                   "CNVPanelizerFromReadCounts"),
#                 envir = environment())
#
#   myCNVPanelizerResults <- parLapply(clust, 1:length(splits), function(i) {
#     CNVPanelizerFromReadCounts(sampleReadCounts[, splits[[i]]],
#                                referenceReadCounts,
#                                genomicRangesFromBed,
#                                numberOfBootstrapReplicates = numberOfBootstrapReplicates,
#                                normalizationMethod = normalizationMethod,
#                                robust = robust,
#                                backgroundSignificanceLevel = backgroundSignificanceLevel,
#                                outputDir = outputDir)
#   })
#   stopCluster(clust)
#
#   n <- length(myCNVPanelizerResults)
#   res <- NULL
#   for (i in seq(n)) {
#     res <- cbind(res, myCNVPanelizerResults[[i]])
#   }
#
#   myCNVPanelizerResultsCompiled <- myCNVPanelizerResults[[1]]
#
#   if(length(myCNVPanelizerResults) > 1) {
#     for (i in 2:length(myCNVPanelizerResults)) {
#       myCNVPanelizerResultsCompiled$sampleReadCounts <- cbind(myCNVPanelizerResultsCompiled$sampleReadCounts, myCNVPanelizerResults[[i]]$sampleReadCounts)
#       myCNVPanelizerResultsCompiled$bootList <- append(myCNVPanelizerResultsCompiled$bootList, myCNVPanelizerResults[[i]]$bootList)
#       myCNVPanelizerResultsCompiled$backgroundNoise <- append(myCNVPanelizerResultsCompiled$backgroundNoise, myCNVPanelizerResults[[i]]$backgroundNoise)
#       myCNVPanelizerResultsCompiled$reportTables <- append(myCNVPanelizerResultsCompiled$reportTables, myCNVPanelizerResults[[i]]$reportTables)
#     }
#   }
#   return(myCNVPanelizerResultsCompiled)
# }

CNVPanelizerFromReadCounts <- function(sampleReadCounts,
                                       referenceReadCounts,
                                       genomicRangesFromBed,
                                       numberOfBootstrapReplicates = 10000,
                                       normalizationMethod = "tmm",
                                       robust = TRUE,
                                       backgroundSignificanceLevel = 0.05,
                                       outputDir = file.path(getwd(), "CNVPanelizer")) {

  #  VerifiyIfOutputDirectoryExistsOrIsNotEmpty(outputDir)


  metadataFromGenomicRanges <- elementMetadata(genomicRangesFromBed)
  geneNames = metadataFromGenomicRanges["geneNames"][, 1]

  #
  # #   metadataFromGenomicRanges <- elementMetadata(genomicRangesFromBed)
  # #   geneNames = metadataFromGenomicRanges["geneNames"][, 1]
  # #   ampliconNames = metadataFromGenomicRanges["ampliconNames"][, 1]
  #
  #   sampleReadCounts <- ReadCountsFromBam(bamFilenames = sampleBamFilepaths,
  #                                         sampleNames = sampleBamFilepaths,
  #                                         gr = genomicRangesFromBed,
  #                                         removeDup = FALSE)
  #
  #   sampleReadCounts <- ReadCountsFromBam(bamFilenames = sampleBamFilepaths,
  #                                         sampleNames = sampleBamFilepaths,
  #                                         gr = genomicRangesFromBed,
  #                                         removeDup = FALSE)

  normalizedReadCounts <- CombinedNormalizedCounts(sampleReadCounts,
                                                   referenceReadCounts,
                                                   method = normalizationMethod
                                                   #                                                   ,                                                   ampliconNames = ampliconNames
  )

  # After normalization data sets need to be splitted again to perform bootstrap
  samplesNormalizedReadCounts = normalizedReadCounts["samples"][[1]]
  referenceNormalizedReadCounts = normalizedReadCounts["reference"][[1]]

  bootList <- BootList(geneNames,
                       samplesNormalizedReadCounts,
                       referenceNormalizedReadCounts,
                       replicates = numberOfBootstrapReplicates)

  # Estimate the background noise left after normalization
  backgroundNoise <- Background(geneNames,
                                samplesNormalizedReadCounts,
                                referenceNormalizedReadCounts,
                                bootList,
                                replicates = numberOfBootstrapReplicates,
                                significanceLevel = backgroundSignificanceLevel,
                                robust = robust)

  # Build report tables
  reportTables <- ReportTables(geneNames,
                               samplesNormalizedReadCounts,
                               referenceNormalizedReadCounts,
                               bootList,
                               backgroundNoise)

  plots = PlotBootstrapDistributions(bootList,
                                     reportTables,
                                     sampleNames = names(reportTables),
                                     outputFolder = file.path(outputDir, "plots"),
                                     save = TRUE)

  WriteListToXLSX(reportTables,
                  multipleFiles = TRUE,
                  outputFolder = file.path(outputDir, "xlsx"))

  # WriteListToXLSX(reportTables, outputFolder =  "", filepath = "list.xlsx") {
  #   write.xlsx(listOfDataFrames, , rowNames = TRUE)
  # }
  #
  #
  #
  #
  # WriteListToXLSX(tmp1@reportTables, filepath = file.path(outputDir, "reportTables.xlsx"))
  # WriteListToXLSX(tmp1@reportTables, filepath = file.path(getwd(), "CNVPanelizer", "reportTables.xlsx"))

  # CNVPanelizerResults <- setClass("CNVPanelizerResults",
  #                                 #                                  setClass("CNVPanelizerResults",
  #
  #                                 # use this function to extract the slot names of the slotNames(results_CNVPanelizer[[1]])
  #                                 slots = c(sampleReadCounts="matrix",
  #                                           referenceReadColunts="matrix",
  #                                           genomicRangesFromBed="GRanges",
  #                                           bootList="list",
  #                                           backgroundNoise="list",
  #                                           plots = "list",
  #                                           reportTables="list"))
  #
  # myCNVPanelizerResults <- CNVPanelizerResults(sampleReadCounts = sampleReadCounts,
  #                                              referenceReadColunts = referenceReadCounts,
  #                                              genomicRangesFromBed = genomicRangesFromBed,
  #                                              bootList = bootList,
  #                                              backgroundNoise = backgroundNoise,
  #                                              plots = plots,
  #                                              reportTables = reportTables)

  myCNVPanelizerResults <- list(sampleReadCounts = sampleReadCounts,
                               referenceReadCounts = referenceReadCounts,
                               genomicRangesFromBed = genomicRangesFromBed,
                               bootList = bootList,
                               backgroundNoise = backgroundNoise,
                               plots = plots,
                               reportTables = reportTables)



  return(myCNVPanelizerResults)
}

CNVPanelizer <- function(sampleBamFilepaths,
                         referenceBamFilepaths,
                         bedFilepath,
                         amplColumnNumber = 6,
                         minimumMappingQuality = 20,
                         numberOfBootstrapReplicates = 10000,
                         removePcrDuplicates = TRUE,
                         #                         analysisMode = "gene",   # analysisMode can be "gene" or "amplicon"
                         robust = TRUE,
                         backgroundSignificanceLevel = 0.05,
                         outputDir = file.path(getwd(), "CNVPanelizer")) {

  #  VerifiyIfOutputDirectoryExistsOrIsNotEmpty(outputDir)

  genomicRangesFromBed <- BedToGenomicRanges(bedFilepath,
                                             ampliconColumn = amplColumnNumber,
                                             split = "_")

  #  metadataFromGenomicRanges <- elementMetadata(genomicRangesFromBed)
  #  geneNames = metadataFromGenomicRanges["geneNames"][, 1]
  #  ampliconNames = metadataFromGenomicRanges["ampliconNames"][, 1]

  sampleReadCounts <- ReadCountsFromBam(bamFilenames = sampleBamFilepaths,
                                        sampleNames = sampleBamFilepaths,
                                        gr = genomicRangesFromBed,
                                        minimumMappingQuality = minimumMappingQuality,
                                        removeDup = removePcrDuplicates)

  referenceReadCounts <- ReadCountsFromBam(bamFilenames = referenceBamFilepaths,
                                           sampleNames = referenceBamFilepaths,
                                           gr = genomicRangesFromBed,
                                           minimumMappingQuality = minimumMappingQuality,
                                           removeDup = removePcrDuplicates)

  results <- CNVPanelizerFromReadCounts(sampleReadCounts = sampleReadCounts,
                                        referenceReadCounts = referenceReadCounts,
                                        genomicRangesFromBed = genomicRangesFromBed,
                                        #                                        bedFilepath = bedFilepath,
                                        #                                        amplColumnNumber = amplColumnNumber,
                                        #                                       minimumMappingQuality = minimumMappingQuality,
                                        numberOfBootstrapReplicates = numberOfBootstrapReplicates,
                                        #                                       removePcrDuplicates = removePcrDuplicates,
                                        #                                       analysisMode = analysisMode,
                                        robust = robust,
                                        backgroundSignificanceLevel = backgroundSignificanceLevel,
                                        outputDir = outputDir)
  return(results)
}


Strict <- function(cnvPanelizerResults,
                   amplificationMinimumThreshold = 2,
                   deletionMaximumThreshold = 0.6,
                   minimumNumberOfAmplicons = 2
                   #                   , fullConsensusAmongAmplicons = TRUE
) {
#  myCNVPanelizerTableResults <- CollectColumnFromAllReportTables(cnvPanelizerResults@reportTables, "ReliableStatus")
  myCNVPanelizerTableResults <- CollectColumnFromAllReportTables(cnvPanelizerResults$reportTables, "ReliableStatus")
  myCNVPanelizerTableResults <- t(myCNVPanelizerTableResults)
  myCNVPanelizerTableResults <- revalueDF(myCNVPanelizerTableResults, c("Normal"=""))

  interestingColumn <- "MeanBoot"
#  meanBootResults <- CollectColumnFromAllReportTables(cnvPanelizerResults@reportTables, interestingColumn)
  meanBootResults <- CollectColumnFromAllReportTables(cnvPanelizerResults$reportTables, interestingColumn)
  class(meanBootResults) <- "numeric"
  meanBootResults <- t(meanBootResults)

#  consensusResults <- ConsensusCheck(cnvPanelizerResults@genomicRangesFromBed, cnvPanelizerResults@sampleReadCounts, cnvPanelizerResults@referenceReadCounts)
  consensusResults <- ConsensusCheck(cnvPanelizerResults$genomicRangesFromBed, cnvPanelizerResults$sampleReadCounts, cnvPanelizerResults$referenceReadCounts)
  consensusResults <- t(consensusResults)

  test_that("colnames from CNVPanelizer results and Consensus match", {
    expect_equal(colnames(consensusResults), colnames(myCNVPanelizerTableResults))
  })

  test_that("rownames from CNVPanelizer results and Consensus match", {
    expect_equal(rownames(consensusResults), rownames(myCNVPanelizerTableResults))
  })

  cNVPanelizerResultsConsensus <- myCNVPanelizerTableResults
  for (rowname in rownames(myCNVPanelizerTableResults)) {
    for (colname in colnames(myCNVPanelizerTableResults)) {
      #        print(paste(rowname, colname))
      if (!consensusResults[rowname, colname]) {
        cNVPanelizerResultsConsensus[rowname, colname] <- ""
        print(paste("Setting ", myCNVPanelizerTableResults[rowname, colname], "to", cNVPanelizerResultsConsensus[rowname, colname]))
      }
    }
  }

  #  table(unlist(myCNVPanelizerTableResults))
  #  table(unlist(cNVPanelizerResultsConsensus))

  test_that("colnames from CNVPanelizer results and meanBootResults match", {
    expect_equal(colnames(meanBootResults), colnames(cNVPanelizerResultsConsensus))
  })

  test_that("rownames from CNVPanelizer results and meanBootResults match", {
    expect_equal(rownames(meanBootResults), rownames(cNVPanelizerResultsConsensus))
  })

  cNVPanelizerResultsConsensusAndThreshold <- cNVPanelizerResultsConsensus
  for (rowname in rownames(cNVPanelizerResultsConsensusAndThreshold)) {
    for (colname in colnames(cNVPanelizerResultsConsensusAndThreshold)) {
      #      print(cNVPanelizerResultsConsensus[rowname, colname])}}
      if ((cNVPanelizerResultsConsensusAndThreshold[rowname, colname] == "Amplification")
          & ((meanBootResults[rowname, colname] <= amplificationMinimumThreshold))
          #          & ((meanBootResults[rowname, colname] <= 3))
      ) {
        cNVPanelizerResultsConsensusAndThreshold[rowname, colname] <- ""
        print(paste("Setting ", cNVPanelizerResultsConsensus[rowname, colname], "to", cNVPanelizerResultsConsensusAndThreshold[rowname, colname], "because", meanBootResults[rowname, colname]))
      }

      if ((cNVPanelizerResultsConsensusAndThreshold[rowname, colname] == "Deletion")
          & ((meanBootResults[rowname, colname] >= deletionMaximumThreshold))
          #& ((meanBootResults[rowname, colname] >= 0.7))
      ) {
        cNVPanelizerResultsConsensusAndThreshold[rowname, colname] <- ""
        print(paste("Setting ", cNVPanelizerResultsConsensus[rowname, colname], "to", cNVPanelizerResultsConsensusAndThreshold[rowname, colname]))
      }
    }
  }

  #  table(unlist(cNVPanelizerResultsConsensus))
  #  table(unlist(cNVPanelizerResultsConsensusAndThreshold))

  cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons <- cNVPanelizerResultsConsensusAndThreshold
#  genesWithMinimumNumberOfAmplicons <- HaveMininumNumberOfAmplicons(cnvPanelizerResults@genomicRangesFromBed, minimumNumberOfAmplicons)
  genesWithMinimumNumberOfAmplicons <- HaveMininumNumberOfAmplicons(cnvPanelizerResults$genomicRangesFromBed, minimumNumberOfAmplicons)
  for (rowname in rownames(cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons)) {
    for (colname in names(genesWithMinimumNumberOfAmplicons)) {
      if ((cNVPanelizerResultsConsensusAndThreshold[rowname, colname] == "Deletion") &
          !genesWithMinimumNumberOfAmplicons[colname]) {
        cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons[rowname, colname] <- ""
        print(paste("Setting ", cNVPanelizerResultsConsensusAndThreshold[rowname, colname], "to", cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons[rowname, colname]))
      }
    }
  }

  #  table(unlist(cNVPanelizerResultsConsensusAndThreshold))
  #  table(unlist(cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons))

  #  write.xlsx(cNVPanelizerResults, file=file.path(outputDir, paste0("cNVPanelizerResults", ".xlsx")), row.names = TRUE)
  #  write.xlsx(cNVPanelizerResultsConsensus, file=file.path(outputDir, paste0("cNVPanelizerResultsConsensus", ".xlsx")), row.names = TRUE)
  #  write.xlsx(cNVPanelizerResultsConsensusAndThreshold, file=file.path(outputDir, paste0("cNVPanelizerResultsConsensusAndThreshold", ".xlsx")), row.names = TRUE)
  #  write.xlsx(cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons, file=file.path(outputDir, paste0("cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons", ".xlsx")), row.names = TRUE)

  return(cNVPanelizerResultsConsensusAndThresholdAndMinimumNumberOfAmplicons)
}


# TODO Could be provided as a measurement of quality for the run
# Overview <- function(
#   #  sampleIdentifier,
#   sampleReadCounts1,
#   referenceReadCounts1,
#   allSamplesReadCounts,
#   genomicRangesFromBed,
#   #  normalizationMethod = "tss",
#   normalizationMethod = "tmm",
#   outputDir = file.path(getwd(), "CNVPanelizer", "overviewKNNReference30")) {
#
#   #  sampleIdentifier <- "I:/Run Data/Run S5-56/M1_S5-56_IonXpress_074_rawlib.bam"
#
#   #   Overview(allSamplesReadCounts[, filepaths],
#   #            allSamplesReadCounts[, referenceReadCounts],
#   #            #         allSamplesReadCounts,
#   #            genomicRangesFromBed,
#   #            #                     normalizationMethod = "tss",
#   #            #                     normalizationMethod = "tmm",
#   #            outputDir = file.path(getwd(), "CNVPanelizer", "overview"))
#
#   sampleReadCounts1 <- allSamplesReadCounts[, filepaths]
#   #  referenceReadCounts1 <- allSamplesReadCounts[, referenceReadCounts]
#   referenceReadCounts1 <- allSamplesReadCounts[, referenceFilepathsKmeans]
#
#   allSamplesReadCountsNormalized <- NormalizeCounts(allSamplesReadCounts, method = normalizationMethod)
#   samplesNormalizedReadCounts = allSamplesReadCountsNormalized[, colnames(sampleReadCounts1)]
#   referenceNormalizedReadCounts = allSamplesReadCountsNormalized[, colnames(referenceReadCounts1)]
#
#   runSamples <- samplesNormalizedReadCounts
#   refSamples <- referenceNormalizedReadCounts
#   entireSetOfSamples <- allSamplesReadCountsNormalized
#
#   geneNames <- unique(genomicRangesFromBed$geneNames)
#   ampliconNames <- genomicRangesFromBed$ampliconNames
#   # TODO validate if ampliconNames are the rownames of samples matrices
#   # ampliconNames <- rownames(allSamplesNormalized)
#
#   sampleIdentifiers <- colnames(runSamples)
#   #sampleIdentifiers <- sampleIdentifiers[1]
#   for (sampleIdentifier in sampleIdentifiers) {
#     myDir <- file.path(outputDir, basename(sampleIdentifier))
#     dir.create(myDir, recursive = TRUE)
#     for (geneName in geneNames) {
#       #    geneName <- "BRAF"
#       myAmpliconsIndexes <- grep(geneName, ampliconNames)
#       png(file=file.path(myDir, paste0(geneName,'_amplicons.png')),
#           width=1000,
#           height=(ceiling(length(myAmpliconsIndexes)/2) * 500))
#
#       #    par(mfrow=c(ceiling(length(myAmpliconsIndexes)/2), 2), mar=c(3,3,5,1))
#       par(mfrow=c(ceiling(length(myAmpliconsIndexes)/2), 2), mar=c(3, 3, 5, 1), oma=c(0,0,7,0))
#
#       for (ampliconName in ampliconNames[myAmpliconsIndexes]) {
#         #  ampliconName <- "chr7_140453102_140453221_BRAF_Ex15-NM_004333"
#         tmp <- data.frame(y = c(as.vector(runSamples[ampliconName, ])
#                                 ,as.vector(refSamples[ampliconName, ])
#                                 ,as.vector(entireSetOfSamples[ampliconName, ])
#         ),
#
#         x = factor(
#           c(rep("runSamples", length(as.vector(runSamples[ampliconName, ])))
#             ,rep("refSamples",length(as.vector(refSamples[ampliconName, ])))
#             ,rep("allSamples",length(as.vector(entireSetOfSamples[ampliconName, ])))
#           ), levels = c("runSamples", "refSamples", "allSamples"))
#
#         #                         x = c(rep("runSamples", length(as.vector(runSamples[ampliconName, ])))
#         #                               ,rep("refSamples",length(as.vector(refSamples[ampliconName, ])))
#         #                               )
#         )
#
#         myValues <- runSamples[ampliconName, sampleIdentifier]
#
#         boxplot(y~x,
#                 data=tmp,
#                 ylab="Read Counts",
#                 xlab="Sample set",
#                 names = c("Run Samples", "Reference", "All Samples"),
#                 #              main = ampliconName)
#                 main = paste(myValues, ampliconName))
#         stripchart(y ~ x,
#                    vertical = TRUE,
#                    data = tmp,
#                    method = "jitter",
#                    add = TRUE,
#                    #                pch = 20,
#                    pch = 20,
#                    col = 'blue')
#
#
#         #    stripchart(list(0.5, 1), vertical=TRUE, add=TRUE, method="stack", col='red', pch="*")
#         stripchart(myValues, vertical=TRUE, add=TRUE, method="stack", col='red', pch=16, cex=2)
#       }
#       title(main = paste(geneName, "-", sampleIdentifier), font.main = 4, outer=TRUE)
#       dev.off()
#     }
#   }
# }

## TODO CHECK THIS
#
#
# #rowIndexesToMerge <- IndexGenesPositions(lungGenomicRanges$geneNames)
# #GeneMultilinePlot(normalSamples, rowIndexesToMerge, normalizationMethod = "tss")
# GeneMultilinePlot <- function(dataSamples, rowIndexesToMerge, normalizationMethod =  "tmm") {
#   dataSamples <- round(NormalizeCounts(dataSamples))
#   result <- NULL
#   for (i in 1:ncol(dataSamples)) {
#     result <- cbind(result, sapply( rowIndexesToMerge, function(x) {mean(dataSamples[x, i])} ))
#   }
#   colnames(result) <- colnames(dataSamples)
#   dataSamples <- result
#   #  dataSamples <- cbind(geneReadCounts = seq_along(rownames(dataSamples)), dataSamples)
#
#   melted <- melt(dataSamples)
#
#   ggplot(data=melted, aes(x=Var1, y=value, group=Var2, colour=Var2, linetype = Var2)) +
#     #geom_line() +
#     #    scale_colour_manual(values = c('pink','orange','white', 'red')) +
#     theme(legend.position="left") +
#     geom_point(aes(color=Var2)) +
#     geom_line(size=1) +
#     theme(axis.text.x = element_text(angle = 90, hjust = 1))
# }
#

# #Example for testing MultiLinePlot
# dataSamples <- read.table(text = "CN CNVPanelizer panelcn.mops ioncopy
#                  CN1 0.1 0.4 0.4
#                  CN2 0.2 0.1 0.3
#                  CN3 0.1 0.2 0.1
#                  CN4 0.1 0.4 0.4
#                  CN5 0.2 0.1 0.3
#                  CN6 0.1 0.3 0.5
#                  CN7 0.6 0.6 0.6
#                  CN8 0.6 0.3 0.2
#                  CN9 0.6 0.8 0.8
#                  CN10 1 0.9 0.8", header=TRUE)

# TODO reuse this code at GeneMultilinePlot ...
# MultiLinePlot <- function(dataSamples, idColumnIdentifier = NULL, title = "", legendPosition = "left", xLabel = "", yLabel = "", smooth = FALSE) {
#   rs <- dataSamples
#   if (is.null(idColumnIdentifier)) {
#     melted = melt(rs)
#   } else {
#     melted = melt(rs, id.vars=idColumnIdentifier)
#   }
#   melted$CN <- as.character(melted$CN)
#   melted$CN <-factor(melted$CN, levels = unique(melted$CN))
#   myPlot <- ggplot(data=melted, aes(x=CN, y=value, group=variable, colour=variable, linetype = variable)) +
#     scale_colour_discrete(name = "Method") + #  scale_colour_manual(values = c('pink','orange','white', 'red')) +
#     ggtitle(title) +
#     geom_point(aes(color=variable)) + # plots the dots..
#     theme(text = element_text(size=10),
#           axis.text.x = element_text(angle=90, hjust=1)) +
#     theme(panel.border = element_blank(),
#           panel.grid.major = element_blank(),
#           panel.grid.minor = element_blank(),
#           axis.line = element_line(colour = "black")) +
#     theme_bw() + # sets the background colour
#     theme(plot.title = element_text(hjust = 0.5)) +
#     theme(legend.position=legendPosition) + xlab(xLabel) + ylab(yLabel)
#   if(smooth) {
#     #      myPlot = myPlot + geom_smooth(method = "loess", se = FALSE, span = 1.5, show.legend = FALSE)
#     myPlot = myPlot + geom_smooth(method = "auto", se = FALSE, span = 1.5, show.legend = FALSE)
#   } else {
#     myPlot = myPlot + geom_line(size = 1, show.legend = FALSE)    # links the dots with lines..
#   }
#   return(myPlot)
# }


SampleReadCountsPools <- function(ampliconsReadCounts, numberOfPools) {
  listOfPools <- list()
  for (i in 1:numberOfPools) {
    listOfPools[[i]] <- ampliconsReadCounts[seq(from = i, to = length(ampliconsReadCounts), by = numberOfPools)]
  }
  return(listOfPools)
}


SampleReadCountsByPool <-function(ampliconReadCounts, pools) {
  # existingPools <- as.numeric(names(table(pools)))
  # if (any(is.na(existingPools)) {
  #   stop("Pools are not defined as 1, 2, ...")
  # }
  existingPools <- names(table(pools))
  if (length(pools) != length(ampliconReadCounts)) {
    stop("Length  of pools is different than ampliconReadCounts")
  }
  listOfPools <- list()
  for (i in existingPools) {
    listOfPools[[i]] <- (ampliconReadCounts[which(pools==i)])
  }
  return(listOfPools)
}


PoolsPlots <- function(sampleReadCountsByPool) {
  plotValues <- as.integer(unlist(lapply(sampleReadCountsByPool, mean)))
  barplot(plotValues, ylab = "Read Counts mean", xlab="Pools", names.arg = names(sampleReadCountsByPool), col = "dodgerblue4",
          #          ylim=c(0, trunc(max(plotValues)))
  )
}








################################################################################
# KaryotypeAberrationPlot
################################################################################
# KaryotypeAberrationPlot(bedFilepath = bedFilepath,
#                         ampliconColumn = 4,
#                         filepath = "TESTANDO_APAGAR_KARYOTYPE.png")

KaryotypeAberrationPlot <- function(bedFilepath,
                                    ampliconColumn,
                                    status = sample(c("red", "green", "blue"), 30, replace = TRUE),
                                    filepath = "KaryotypeAberrationPlot.png",
                                    width = 1200,
                                    height = 800) {
  #if (!requireNamespace("BiocManager", quietly=TRUE))
      #install.packages("BiocManager")
  #BiocManager::install("chromPlot")
  #library(chromPlot)
  #library(stringr)
  hg_gap <- NULL # not sure if this is good solution to avoid 'NOTE'
  data(hg_gap)
  getGeneRegionsStatusSummary <- GetGeneRegionsStatusSummary(bedFilepath = bedFilepath,
                                                             ampliconColumn = ampliconColumn,
                                                             status = sample(c("red", "green", "blue"), 30, replace=TRUE))
  geneRegionsBands <- GetGeneRegionsBand(getGeneRegionsStatusSummary)
  png(filepath, width = width, height = height)
  #  chromPlot(gaps=hg_gap)
  #  chromPlot(gaps=hg_gap, stat=geneRegions, statCol="Value")
  #  chromPlot(gaps=hg_gap, stat=geneRegions, statCol="Value", chr=unique(geneRegions$Chrom))
  #  chromPlot(gaps=hg_gap, bands=hg_cytoBandIdeo, stat=geneRegions, statCol="Value", chr=unique(geneRegions$Chrom), statName="Value", noHist=TRUE, figCols=5, cex=0.7, statTyp="n", chrSide=c(1,1,1,1,1,1,-1,1))

  # TODO TEMPORARARIO PORQUE O chromPlot nao pode ser instalado em versao 3.2.0 ...
  # karyotypeAberrationPlot <- chromPlot(gaps = hg_gap,
  #                                      bands = geneRegionsBands,
  #                                      chr = unique(getGeneRegionsStatusSummary$Chrom),
  #                                      stat = getGeneRegionsStatusSummary,
  #                                      statCol = "Value",
  #                                      statName = "Value",
  #                                      noHist = TRUE,
  #                                      figCols = 7,
  #                                      cex = 0.7,
  #                                      statTyp = "n",
  #                                      chrSide = c(1,1,1,1,1,1,-1,1))
  karyotypeAberrationPlot <- NA
  dev.off()
  return(karyotypeAberrationPlot)
}


# It Merges all regions related to the same gene and assigns it a color that relates to the aberratiosn status
# Status should be a vector of colors representing the amplifications (red for amplification, green for deletion and blue for normal)
GetGeneRegionsStatusSummary <- function(bedFilepath, ampliconColumn = 4, split = "_", status = c("blue")) {
  a <- read.csv(bedFilepath, sep="\t", skip=1, header = FALSE)
  a[, 1] <- str_replace(a$V1, "chr","")
  a[, ampliconColumn] <- vapply(strsplit(as.vector(a[, ampliconColumn]), "_", fixed = TRUE), "[", "", 1)
  a <- a[,1:4]
  colnames(a) <- c("Chrom", "Start", "End", "ID"
                   #                   , "Colors"
  )
  geneNames <- unique(a$ID)
  geneRegions <- NULL
  for (geneName in geneNames) {
    genomicRegions <- a[a$ID==geneName, ]
    chromossome <- unique(genomicRegions$Chrom)
    start <- min(genomicRegions$Start)
    end <- max(genomicRegions$End)
    gene <- unique(genomicRegions$ID)
    geneRegions <- rbind(geneRegions, c(chromossome, start, end, gene))
  }
  colnames(geneRegions) <- colnames(a)
  geneRegions <- as.data.frame(geneRegions,
                               #                               as.is = TRUE)
                               stringsAsFactors = FALSE)
  geneRegions <- cbind(geneRegions, Colors = status)

  #  geneRegions$Chrom <- as.numeric(geneRegions$Chrom)
  geneRegions$Start <- as.numeric(geneRegions$Start)
  geneRegions$End <- as.numeric(geneRegions$End)
  geneRegions$ID <- as.factor(geneRegions$ID)
  return(geneRegions)
}


# Returns objected needed for plotting.. it changes some column names and makes the regions larger to be able to be visible in the plot..
#GetGeneRegionsBand <- function(geneRegionsStatusSummary) {
GetGeneRegionsBand <- function(geneRegions) {
  geneRegionsBands <- geneRegions
  colnames(geneRegionsBands) <- c("Chrom", "Start", "End", "Name", "Colors")
  geneRegionsBands$Chrom <- paste0("chr", geneRegionsBands$Chrom)
  geneRegionsBands$Name <- as.character(geneRegionsBands$Name)
  geneRegionsBands$Colors <- as.character(geneRegionsBands$Colors)
  # since the regions are too small to be visualized.. around 1 milion times smaller than the size of chromossome..
  geneRegionsBands$End <- geneRegionsBands$End + 1E6
  return(geneRegionsBands)
}






# TODO
# Could take in consideration the size of regions and
# the whole coverage of the gene (number of amplicons * seqLength)
# > 50 bp
# panelStatistics(genomicRangesFromBed) {
#   # hist(width(ranges(genomicRangesFromBed)))
#   # genomicRangesFromBed$geneNames
#
# }



#   boxplotGeneLevel <- function(allSamplesNormalized,
#                                cohortFilepaths,
#                                referenceFilepathsKNN,
#                                referenceFilepathsPercentil,
#                                geneName,
#                                outputDir = file.path(getwd(),"batchAnalysisNoMininumMappingQuality2")) {
#     #boxplotPerAmplicon <- function(ampliconIndex) {
#     #geneName <- "BRAF"
#     #geneName <- "APC"
#     #  geneName <- "ERBB4"
#
# #    setwd(outputDir)
#
#     #  cohortFilepaths <- filepaths
#     martina <- allSamplesNormalized[, (colnames(allSamplesNormalized) %in% cohortFilepaths)]
#
#     nonMartina <- allSamplesNormalized[, !(colnames(allSamplesNormalized) %in% cohortFilepaths)]
#
#
#     for (geneName in geneNames) {
#
#
#
#       ampliconNames <- rownames(allSamplesNormalized)
#       myAmpliconsIndexes <- grep(geneName, ampliconNames)
#
#       #  pdf(file=paste0(geneName,'_plot.pdf'), width=10, height= ceiling(length(myAmpliconsIndexes)/2) * 5)
#       png(file=paste0(geneName,'_plot.png'), width=10, height= ceiling(length(myAmpliconsIndexes)/2) * 5)
#
#       par(mfrow=c(ceiling(length(myAmpliconsIndexes)/2), 2), mar=c(3,3,5,1))
#
#       for (ampliconName in ampliconNames[myAmpliconsIndexes]) {
#         #  ampliconName <- "chr7_140453102_140453221_BRAF_Ex15-NM_004333"
#         tmp <- data.frame(y = c(as.vector(nonMartina[ampliconName, ]),
#                                 as.vector(nonMartina[ampliconName, referenceFilepathsKNN]),
#                                 as.vector(nonMartina[ampliconName, referenceFilepathsPercentil]),
#                                 as.vector(martina[ampliconName, ])),
#                           x = c(rep("nonMartina",length(as.vector(nonMartina[ampliconName, ]))),
#                                 rep("referenceKNN", length(as.vector(nonMartina[ampliconName, referenceFilepathsKNN]))),
#                                 rep("referencePercentil", length(as.vector(nonMartina[ampliconName, referenceFilepathsPercentil]))),
#                                 rep("martina", length(as.vector(martina[ampliconName, ])))))
#
#         boxplot(y~x, data=tmp, ylab="Read Counts", xlab="Sample set", main = ampliconName)
#         stripchart(y ~ x,
#                    vertical = TRUE,
#                    data = tmp,
#                    method = "jitter",
#                    add = TRUE,
#                    pch = 20,
#                    col = 'blue')
#       }
#       dev.off()
#     }
#   }





# i <- 0
# for (x in 1:100) {
#   i <- i+1
#   if (i %% 10 == 0) {
#     Sys.sleep(1)
#     cat(".")
#   }
# }
#