In adamsardar/metaDEGth: Combining differential expression P-values in a statistically consistent fashion whilst respecting background signal/noise trends

erlotinib_adrHCC827_diffexDT

susceptibleCellLine_BetaUniformModel <- fitBetaUniformMixtureDistribution(erlotinib_adrHCC827_diffexDT$HCC827_pValue, 
                                                                          nStarts = 20)
fittedBetaShape <- susceptibleCellLine_BetaUniformModel@a
uniformProportion <- susceptibleCellLine_BetaUniformModel@lambda

fitBetaUniformMixtureDistribution(rbetauniform(100000,fittedBetaShape,nonUniformProportion))

plot.BetaUniformModel(susceptibleCellLine_BetaUniformModel, outputParameters = F, outputFormula = F)

qplot(rbetauniform(10000, a = susceptibleCellLine_BetaUniformModel@a, lambda = susceptibleCellLine_BetaUniformModel@lambda))



nullModelTests <- list()

for(i in 1:1E4){

  randomSetSize <- rpois(1,10)

  random_draws <- map(1:1E4, ~{rbetauniform(randomSetSize, fittedBetaShape, uniformProportion)})
  # random_draws <- map(1:1000, ~{runif(n = randomSetSize)}); uniformProportion <- 1
  #random_draws <- map(1:1000, ~{rbeta(n = randomSetSize, shape1 = fittedBetaShape, shape2 = 1)}) ; uniformProportion <- 0

  testStats <- map_dbl(random_draws, ~ sum(-log(.x)))

  empiricalGTE <- length(which(testStats[1] <= testStats[-1]))/length(testStats[-1])

  fisherPval <- pgamma(testStats[1], shape = randomSetSize, scale = 1, lower.tail = FALSE) 

  mu <- randomSetSize*(uniformProportion +  (1-uniformProportion)*scaleParam)
  sigmaSq <-  randomSetSize*(uniformProportion +  (1-uniformProportion)*scaleParam^2)

  combinedPval <- pgamma(testStats[1], shape = randomSetSize*uniformProportion, scale = 1, lower.tail = FALSE) + 
    pgamma(testStats[1], shape = randomSetSize*(1-uniformProportion), rate = fittedBetaShape, lower.tail = FALSE)

  # hyperexponentialSumTransitionMatrix <- constructHyperexponentialSumTransitionMatrix(randomSetSize, uniformProportion, fittedBetaShape)
  # 
  # initialProb <- Matrix(0,ncol = ncol(hyperexponentialSumTransitionMatrix), nrow = 1)
  # initialProb[1:2] <- c(uniformProportion, 1-uniformProportion)

  # CDF of phase rate distribution os 1- alpha * exp(x*S)*oneVec
 # phaseRateP <- sum(initialProb %*% Matrix::expm(testStats[1]*hyperexponentialSumTransitionMatrix))

  phaseRateP <- betaUniformPvalueSumTest(random_draws[[1]], susceptibleCellLine_BetaUniformModel)
  # Welch-Satterthwaite Approximation for a sum of exponentials, some with different paramenters
  approxPval <- pgamma(testStats[1], shape = (mu^2)/sigmaSq, scale = sigmaSq/mu, lower.tail = F)

  nullModelTests %<>% c(list(data.table(empiricalGTE, fisherPval, combinedPval, approxPval, phaseRateP, size = randomSetSize)))
}


microbenchmark::microbenchmark(

    combinedPval <- pgamma(testStats[1], shape = randomSetSize*uniformProportion, scale = 1, lower.tail = FALSE) + 
    pgamma(testStats[1], shape = randomSetSize*(1-uniformProportion), rate = fittedBetaShape, lower.tail = FALSE),

  hyperexponentialSumTransitionMatrix <- constructHyperexponentialSumTransitionMatrix(randomSetSize, uniformProportion, fittedBetaShape),

  phaseRateP <- pphtype(testStats[1],
                        prob = c(c(uniformProportion, 1-uniformProportion), rep(0,2*(randomSetSize-1))),
                        rates = as.matrix(hyperexponentialSumTransitionMatrix),
                        lower.tail = FALSE),

  # Welch-Satterthwaite Approximation for a sum of exponentials, some with different paramenters
  approxPval <- pgamma(testStats[1], shape = (mu^2)/sigmaSq, scale = sigmaSq/mu, lower.tail = F)
)

ggplot(data = data.frame(x = c(0, 100)), aes(x)) +

  stat_function(fun = dgamma, args = list(rate = fittedBetaShape, shape = randomSetSize), colour = "red") +
  stat_function(fun = dgamma, args = list(rate = 1, shape = randomSetSize), colour = "black") +
  stat_function(fun = dgamma, args = list(scale = sigmaSq/mu, shape = (mu^2)/sigmaSq), colour = "pink") +

  geom_histogram(data = data.table(x = testStats), aes( y = ..density..), binwidth = 1, fill = "grey", color = "black") +

  stat_function(fun = actuar::dphtype, args = list(prob = c(c(uniformProportion, 1-uniformProportion), rep(0,2*(randomSetSize-1))),
                                           rates = as.matrix(hyperexponentialSumTransitionMatrix)), colour = "green", size = 1) +

  ylab("") +
  scale_y_continuous(breaks = NULL) +
  geom_vline(xintercept = testStats[1]) 


empiricalVsPred <- rbindlist(nullModelTests)

empiricalVsPred[empiricalGTE <= 1E-2]



empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE - combinedPval))] + 
  scale_x_log10(breaks = 10^(-seq(0,6))) +
  scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) +
  theme_bw()

empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE - approxPval))] + 
  scale_x_log10(breaks = 10^(-seq(0,6))) +
  scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) +
  theme_bw()

empiricalVsPred[,qplot(x = combinedPval, y = (empiricalGTE -  phaseRateP))] + 
  scale_x_log10(breaks = 10^(-seq(0,6))) +
  scale_y_continuous(breaks = seq(-0.13,0.1,0.01)) +
  theme_bw()


empiricalVsPred[,qplot(x = size, y = (empiricalGTE - phaseRateP))]

empiricalVsPred[empiricalGTE <= 1E-2] %>%
  ggplot(aes(x = empiricalGTE)) +
  geom_point(aes(y = fisherPval), colour = "blue") +
  geom_point(aes(y = combinedPval), colour = "red") +
  geom_point(aes(y = approxPval), colour = "green") +
  geom_point(aes(y = phaseRateP), colour = "pink") +
  geom_abline() +
  theme_bw() +
  scale_x_log10() +
  scale_y_log10()

adamsardar/metaDEGth documentation built on Sept. 28, 2022, 10:12 p.m.