#Calculation of quantile normalization factors necessary for calcNormFactorsGSD function. Same for both train and test sets.

calcFactorQuantileGSD =

  function (data, lib.size, p = 0.75) 

  {

    y <- t(t(data)/lib.size)

    f <- apply(y, 2, function(x) quantile(x, p = p))

  }

#Calculation of weighted normalization factors necessary for calcNormFactorsGSD function.

calcFactorWeightedGSD =

  function (obs, ref, libsize.obs = NULL, libsize.ref = NULL, logratioTrim = 0.3, 

            sumTrim = 0.05, doWeighting = TRUE, Acutoff = -1e+10) 

  {

    if (all(obs == ref)) 

      return(1)

    obs <- as.numeric(obs)

    ref <- as.numeric(ref)

    if (is.null(libsize.obs)) 

      nO <- sum(obs)

    else nO <- libsize.obs

    if (is.null(libsize.ref)) 

      nR <- sum(ref)

    else nR <- libsize.ref

    logR <- log2((obs/nO)/(ref/nR))

    absE <- (log2(obs/nO) + log2(ref/nR))/2

    v <- (nO - obs)/nO/obs + (nR - ref)/nR/ref

    fin <- is.finite(logR) & is.finite(absE) & (absE > Acutoff)

    logR <- logR[fin]

    absE <- absE[fin]

    v <- v[fin]

    n <- sum(fin)

    loL <- floor(n * logratioTrim) + 1

    hiL <- n + 1 - loL

    loS <- floor(n * sumTrim) + 1

    hiS <- n + 1 - loS

    keep <- (rank(logR) >= loL & rank(logR) <= hiL) & (rank(absE) >= loS & rank(absE) <= hiS)

    if (doWeighting) 

      2^(sum(logR[keep]/v[keep], na.rm = TRUE)/sum(1/v[keep],na.rm = TRUE))

    else 2^(mean(logR[keep], na.rm = TRUE))

  }

#Calculation of deseq normalization factors necessary for calcNormFactorsGSD function.

calcFactorRLEGSD =

  function (data.train, data.test, lib.size, lib.size.test) 

  {

    gm <- exp(rowMeans(log(data.train)))

    f = apply(data.train, 2, function(u) median((u/gm)[gm > 0]))

    f.test = apply(data.test, 2, function(u) median((u/gm)[gm > 0]))

    f = f / lib.size

    f.test = f.test / lib.size.test

    deseqsizefactors = list(f,f.test)

    return(deseqsizefactors)

  }

#Calculation of normalization factors using calcNormFactorsGSD function.

calcNormFactorsGSD =

  function (data.train, data.test, lib.size = NULL, method = c("TMM", "deseq", "none"), refColumn = NULL, logratioTrim = 0.3, sumTrim = 0.05, 

            doWeighting = TRUE, Acutoff = -1e+10, p = 0.75, ...) 

  {

    x <- as.matrix(data.train)

    xtest <- as.matrix(data.test)

    if (any(is.na(x)||is.na(xtest)))

      stop("NAs not permitted")

    if (is.null(lib.size)) 

      lib.size <- colSums(x)

    lib.size.test <- colSums(xtest)

    method <- match.arg(method)

    allzero <- rowSums(x > 0) == 0

    if (any(allzero)) 

      x <- x[!allzero, , drop = FALSE]

      xtest <- xtest[!allzero, , drop = FALSE]

    if (nrow(x) == 0 || ncol(x) == 1) 

      method = "none"

    if(method == "TMM"){

      f75 <- calcFactorQuantileGSD(data = x, lib.size = lib.size, p = 0.75)

      f75.test <- calcFactorQuantileGSD(data = xtest, lib.size = lib.size.test, p = 0.75)

      refColumn <- which.min(abs(f75 - mean(f75)))

      f <- rep(NA, ncol(x))

      f.test <- rep(NA, ncol(xtest))

      for (i in 1:ncol(x)) f[i] <- calcFactorWeightedGSD(obs = x[,i], ref = x[, refColumn], libsize.obs = lib.size[i], 

                                                         libsize.ref = lib.size[refColumn], logratioTrim = logratioTrim, 

                                                         sumTrim = sumTrim, doWeighting = doWeighting, Acutoff = Acutoff)

      for (i in 1:ncol(xtest)) f.test[i] <- calcFactorWeightedGSD(obs = xtest[,i], ref = x[, refColumn], libsize.obs = lib.size.test[i], 

                                                                  libsize.ref = lib.size[refColumn], logratioTrim = logratioTrim, 

                                                                  sumTrim = sumTrim, doWeighting = doWeighting, Acutoff = Acutoff)

      normf = list(f,f.test)

    } 

    else if(method == "deseq"){

      normf = calcFactorRLEGSD(data.train = x, data.test = xtest, lib.size = lib.size, lib.size.test = lib.size.test)#/lib.size 

    }

    else {

      normf = list(rep(1, ncol(x)), rep(1, ncol(xtest)))

    }

    names(normf) = c("train", "test")

    f = as.numeric(normf[[1]]) / (exp(mean(log(normf[[1]]))))

    f.test = as.numeric(normf[[2]]) / (exp(mean(log(normf[[1]]))))

    normf2 = list(f, f.test, lib.size, lib.size.test)

    names(normf2) = c("TrainNormFactor","TestNormFactor","TrainLibSize","TestLibSize")

    return(normf2)

  }

#traindatayı, testdatayı ve trainclassı alacak, method = TMM, RLE, none olacak.

#train ve test icin voom gibi expression ve weights döndürecek.

voomGSD = 

  function(data.train, data.test, group, norm.method = c("TMM", "deseq", "none"), design = NULL, lib.size = NULL, span = 0.5)

  {

    out <- list()

    NormFactors = calcNormFactorsGSD(data.train = data.train, data.test = data.test, method = norm.method)

    TrainNormFactor = NormFactors$TrainNormFactor

    TestNormFactor = NormFactors$TestNormFactor

    TrainLibSize = NormFactors$TrainLibSize

    TestLibSize = NormFactors$TestLibSize

    lib.size.tr = TrainNormFactor * TrainLibSize

    lib.size.ts = TestNormFactor * TestLibSize

    design.tr = model.matrix(~group)

    rownames(design.tr) = colnames(data.train)

    design.ts <- matrix(1, ncol(data.test), 1)

    rownames(design.ts) <- colnames(data.test)

    colnames(design.ts) <- "GrandMean"

    y.tr <- t(log2(t(data.train + 0.5)/(lib.size.tr + 1) * 1e+06))

    y.ts <- t(log2(t(data.test + 0.5)/(lib.size.ts + 1) * 1e+06))

    fit.tr <- lmFit(y.tr, design.tr)

    fit.ts <- lmFit(y.ts, design.ts)

    if (is.null(fit.tr$Amean)) 

      fit$Amean <- rowMeans(y.tr, na.rm = TRUE)

    fit.ts$Amean = fit.tr$Amean

    fit.ts$sigma = fit.tr$sigma

    fit.ts$coefficients = fit.tr$coefficients[,1]

    sx <- fit.tr$Amean + mean(log2(lib.size.tr + 1)) - log2(1e+06)

    sy <- sqrt(fit.tr$sigma)

    l <- lowess(sx, sy, f = span)

    f <- approxfun(l, rule = 2)

    fitted.values.tr <- fit.tr$coefficients %*% t(fit.tr$design)

    fitted.values.ts <- fit.ts$coefficients %*% t(fit.ts$design)

    fitted.cpm.tr <- 2^fitted.values.tr

    fitted.cpm.ts <- 2^fitted.values.ts

    fitted.count.tr <- 1e-06 * t(t(fitted.cpm.tr) * (lib.size.tr + 1))

    fitted.count.ts <- 1e-06 * t(t(fitted.cpm.ts) * (lib.size.ts + 1))

    fitted.logcount.tr <- log2(fitted.count.tr)

    fitted.logcount.ts <- log2(fitted.count.ts)

    w.tr <- 1/f(fitted.logcount.tr)^4

    w.ts <- 1/f(fitted.logcount.ts)^4

    dim(w.tr) <- dim(fitted.logcount.tr)

    dim(w.ts) <- dim(fitted.logcount.ts)

    dimnames(w.tr) = dimnames(y.tr)

    dimnames(w.ts) = dimnames(y.ts)

    out$TrainExp <- y.tr

    out$TestExp <- y.ts

    out$TrainWeights <- w.tr

    out$TestWeights <- w.ts

    new("EList", out)

  }