badge.gamma <- function(data, gene.names, n.a,n.b,t.val){
  ### bayesian differential analysis of gene expression data
  ### works with mixture of log-normal and gamma densityes

   n.genes <- nrow(data)    # defines number of genes

#### Group 2
    # "Analysis of genes that follow Gamma distribution"
        normal.data <- as.matrix(data[,1:n.a])
        tumor.data  <- as.matrix(data[,(n.a+1):(n.a+n.b)])
      normal.group.2 <- normal.data
       tumor.group.2 <- tumor.data
       n.genes.2 <- nrow(normal.group.2)


### compute moment estimators of alpha and beta

      mean.a.hat <- apply(normal.group.2,1,mean)
      mean.b.hat <- apply(tumor.group.2,1,mean)
         alpha.a.hat <- (mean.a.hat)^2/ apply((normal.group.2),1,var)
        alpha.b.hat <- (mean.b.hat)^2/ apply((tumor.group.2),1,var)
 ### prior hyper-parameters

       b.cond1.0 <- 20*n.a; a.cond1.0 <- 1
       b.cond2.0 <- 20*n.b; a.cond2.0 <- 1

 # posterior parameters
     beta.a.pos  <- b.cond1.0+apply(normal.group.2,1,sum);
     alpha.a.pos <- apply( cbind( ncol(normal.group.2)*alpha.a.hat + a.cond1.0,75),1,min)
     beta.b.pos  <-  b.cond2.0+apply(tumor.group.2,1,sum)
     alpha.b.pos <- apply(cbind( ncol(tumor.group.2)* alpha.b.hat + a.cond2.0, 75),1,min)
     theta.gamma <-  (beta.a.pos/n.a)/(beta.b.pos/n.b)*(1+1/(n.b*alpha.b.hat))
 
 # computation with F distribution:

  f.val <- ( beta.b.pos/beta.a.pos) *( n.a+1/alpha.a.hat)/(n.b+1/alpha.b.hat)*(t.val)
    l.ci.f <-  qf(0.025, (2*alpha.b.pos),(2*alpha.a.pos)) * theta.gamma
    u.ci.f <-  qf(0.975, (2*alpha.b.pos),(2*alpha.a.pos)) * theta.gamma
    prob.f <- 1-pf(f.val,(2*alpha.b.pos),(2*alpha.a.pos))
    flag.gamma<- rep(-1,n.genes.2)
 res.gamma <- cbind( gene.names,theta.gamma,prob.f,l.ci.f,u.ci.f,mean.a.hat,alpha.a.hat,
                     mean.b.hat,alpha.b.hat,f.val)
     return(res.gamma)}


 badge.lnorm.dv <- function(data, gene.names, n.a,n.b){
  ### bayesian differential analysis of gene expression data
  ### read functions
   n.genes <- nrow(data)    # defines number of genes
       normal.data <- as.matrix(data[,1:n.a])
        tumor.data  <- as.matrix(data[,(n.a+1):(n.a+n.b)])  
     normal.group.1 <- normal.data
     tumor.group.1  <- tumor.data
     n.norm <- max(ncol(normal.group.1),4)
     n.tumo <- max(ncol(tumor.group.1),4)
     n.genes.1 <- nrow(normal.group.1)
    
     mean.normal  <- apply(  log(normal.group.1),1,mean);
     mean.tumor   <-  apply(  log(tumor.group.1),1,mean);
     var.norm     <-  ( apply(  log(normal.group.1) ,1,var))
     var.tumo     <-  ( apply(  log(tumor.group.1) ,1,var))
     noncentral.t <- 1/2*(var.norm-var.tumo)
     theta.lnorm <- exp(mean.normal -mean.tumor +noncentral.t)
  
  cos.phi <- (var.tumo/n.tumo)/(var.tumo/n.tumo+var.norm/n.norm )
  sin.phi <- 1-cos.phi
  f.1 <- (n.tumo-1)/(n.tumo-3)*cos.phi + (n.norm-1)/(n.norm-3)*sin.phi
     ifelse( n.tumo > 5 , 
          { f.2 <- (n.tumo-1)^2/( (n.tumo-3)^2*(n.tumo-5)) *cos.phi^2 +
                   (n.norm-1)^2/( (n.norm-3)^2*(n.norm-5)) *sin.phi^2;
           b <- 4+f.1^2/f.2}, 
           b <- 4)
  a <- f.1*(b-2)/b
   t.stat.new <- (-noncentral.t - mean.normal +mean.tumor )/
    sqrt(a *(var.norm /(n.norm)+var.tumo /(n.tumo)))
 
  #### now compute prob theta=mean.normal/mean.tumor>1, and 95% credible interval

  prob.t <- 1-pt(t.stat.new , b ) 

  l.ci.t  <- mean.normal-mean.tumor+qt(0.025,b)*sqrt(a*(var.norm/(n.norm)+var.tumo/(n.tumo))) + noncentral.t 
  u.ci.t  <- mean.normal-mean.tumor+qt(0.975,b)*sqrt(a*(var.norm/(n.norm)+var.tumo/(n.tumo))) + noncentral.t 
  res.lnorm <- cbind(gene.names,theta.lnorm,prob.t,exp(l.ci.t),exp(u.ci.t),mean.normal,var.norm,
                      mean.tumor,var.tumo,a,b)

   return(res.lnorm)}

 # function to compute bayes factor
  bayes.factor.ratio <- function(data){
      # prior on tau,mu is tau^(-g)
      # alpha.0 and beta.0 are hyper-parameters of beta prior

       data <- unlist(data) 
       n.col <- length(data)
              
 # posterior parameters
     alpha.hat <-   mean(data)^2/var(data)  
     beta.hat   <-  (alpha.hat+1/n.col)/(mean(data)) 
     logmean    <- mean(log( data ))
     logvar     <- var(log(data))*(n.col-1)/n.col; 
 
    score.gamma <-  (alpha.hat-1)*sum(log(data))-n.col*lgamma(alpha.hat)+ n.col*alpha.hat*log(beta.hat)-
      beta.hat*sum(data)
   

  score.lnorm <-  -(n.col/2)*log(2*pi*logvar)-sum(log(data)) -n.col/2 
      res <- c( score.gamma,score.lnorm)
      return(res)
}


  posterior.weights <- function(data){

    weights <- apply(data,1,bayes.factor.ratio)
    return(weights)}     
 
 ttest.log <- function(data, gene.names, n.a,n.b){
  ### standard t-test
  ### read functions
   n.genes <- nrow(data)    # defines number of genes
        normal.data <- as.matrix(data[,1:n.a])
        tumor.data  <- as.matrix(data[,(n.a+1):(n.a+n.b)])

     normal.group.1 <- normal.data
     tumor.group.1  <- tumor.data

  n.norm <- ncol(normal.group.1)
  n.tumo <- ncol(tumor.group.1)
  n.genes.1 <- nrow(normal.group.1)
    
     mean.normal  <- apply(  log(normal.group.1),1,mean);
     mean.tumor   <-  apply(  log(tumor.group.1),1,mean);
     var.norm     <-  ( apply(  log(normal.group.1) ,1,var))
     var.tumo     <-  ( apply(  log(tumor.group.1) ,1,var))
     t.stat.new <- (mean.normal -mean.tumor )/
          sqrt( var.norm /(n.norm)+var.tumo /(n.tumo))
 
  #### now compute prob theta=mean.normal/mean.tumor>1, and 95% credible interval
   return(t.stat.new)}

 ttest <- function(data, gene.names, n.a,n.b){
  ### standard t-test
  ### read functions
   n.genes <- nrow(data)    # defines number of genes
        normal.data <- as.matrix(data[,1:n.a])
        tumor.data  <- as.matrix(data[,(n.a+1):(n.a+n.b)])

     normal.group.1 <- normal.data
     tumor.group.1  <- tumor.data


#### Group 1: here, first check if log-transformed data have the same variance;
  n.norm <- ncol(normal.group.1)
  n.tumo <- ncol(tumor.group.1)
  n.genes.1 <- nrow(normal.group.1)
    
     mean.normal  <- apply( (normal.group.1),1,mean);
     mean.tumor   <-  apply((tumor.group.1),1,mean);
     var.norm     <-  ( apply( (normal.group.1) ,1,var))
     var.tumo     <-  ( apply( (tumor.group.1) ,1,var))
       t.stat.new <- (mean.normal -mean.tumor )/
    sqrt( var.norm /(n.norm)+var.tumo /(n.tumo))
 
  #### now compute prob theta=mean.normal/mean.tumor>1, and 95% credible interval
    return(t.stat.new)}