"sagat2" <-
function(data,U,D=NULL,nEG=ncol(as.matrix(U)),use.cor=TRUE,df.prior=0)
{
  ## At first I'll assume that incoming data are only variates of the contrast of interest

  # Preliminaries
    U <- U[,1:nEG]
    U2 <- U^2
    U <- as(U,"dgeMatrix")
    U2 <- as(U2,"dgeMatrix")
    if(nrow(data) != nrow(U)) stop("Dimensions of data and U aren't compatible!")
    N <- nrow(data)
    M <- nEG
    m <- ncol(data)

  # Impute missing data if necessary
    if(any(is.na(data))) {
	require(EMV,quietly=TRUE)
        is.row.na <- apply(data,1,function(x) all(is.na(x)))
        data[is.row.na,] <- matrix(colMeans(data,na.rm=TRUE),nrow=sum(is.row.na),ncol=ncol(data),byrow=TRUE)
        data[!is.row.na,] <- knn(data[!is.row.na,],k=10)[[1]]
    }

  # Estimate shared and unique covariance matrices
    Wx1 <- solve(crossprod(U2) - Diagonal(x=rep(1,M)))

    if(use.cor) {

      diagS.vec <- apply(data,1,var)
      Wx0 <- as((t(U) %*% Diagonal(x=sqrt(diagS.vec)) %*% U)^2,"dsyMatrix")
      Wx2 <- rowSums(t(U2))
      R <- as(cor(t(data)),"dsyMatrix")
      Wx3 <- t(U) %*% R; rm(R)
      Wx3 <- sapply(1:nrow(Wx3),function(i) Wx3[i,] %*% U[,i])
      
      WxR <- Wx1 %*% (Wx2 - Wx3); rm(Wx1)
      
      WxR <- as.matrix(WxR)
      WxR[WxR <= 0] <- min(WxR[WxR > 0])
      WxR <- as(WxR,"dgeMatrix")
      
      Wx <- Wx0 %*% WxR; rm(Wx0)

      WeR <- rep(1,N) - (U2 %*% WxR)
      
      WeR <- as.matrix(WeR)
      WeR[WeR <= 0] <- min(WeR[WeR > 0])
      WeR <- as(WeR,"dgeMatrix")
      
      We <- diagS.vec * WeR
    }
    else {

      S <- as(cov(t(data)),"dsyMatrix")
      diagS.vec <- diag(S)
      Wx2 <- t(U2) %*% diagS.vec
      Wx3 <- t(U) %*% S; rm(S)
      Wx3 <- sapply(1:nrow(Wx3),function(i) Wx3[i,] %*% U[,i])

      Wx <- Wx1 %*% (Wx2 - Wx3); rm(Wx1)

      Wx <- as.matrix(Wx)
      Wx[Wx <= 0] <- min(Wx[Wx > 0])
      Wx <- as(Wx,"dgeMatrix")

      We <- diagS.vec - (U2 %*% Wx)

      We <- as.matrix(We)
      We[We <= 0] <- min(We[We > 0])
      We <- as(We,"dgeMatrix")
    }
    if(!is.null(D)) {
      We0 <- sum(D[(nEG+1):length(D)]^2/(length(D)-1)) / (length(D)-nEG)
      We <- (df.prior*We0 + m*We) / (df.prior + m)
    }

  # Return z scores  
    return(rowMeans(data) / sqrt(as.numeric(We/m)))
}