#' Normal mixture model estimation and selection for a series of cluster numbers
#'
#' Perform co-expression and co-abudance analysis of high-throughput
#' sequencing data, with or without data transformation, using a Normal
#' mixture models. The output of \code{NormMixClus} is an S4 object of
#' class \code{RangedSummarizedExperiment}.
#'
#' @param y_profiles (\emph{n} x \emph{q}) matrix of observed profiles for \emph{n}
#' observations and \emph{q} variables
#' @param K Number of clusters (a single value or a sequence of values).
#' @param subset Optional vector providing the indices of a subset of
#' genes that should be used for the co-expression analysis (i.e., row indices
#' of the data matrix \code{y}.
#' @param parallel If \code{FALSE}, no parallelization. If \code{TRUE}, parallel
#' execution using BiocParallel (see next argument \code{BPPARAM}). A note on running
#' in parallel using BiocParallel: it may be advantageous to remove large, unneeded objects
#' from the current R environment before calling the function, as it is possible that R's
#' internal garbage collection will copy these files while running on worker nodes.
#' @param BPPARAM Optional parameter object passed internally to \code{bplapply} when
#' \code{parallel=TRUE}. If not specified, the parameters last registered with \code{register}
#' will be used.
#' @param seed If desired, an integer defining the seed of the random number generator. If
#' \code{NULL}, a random seed is used.
#' @param ... Additional optional parameters to be passed to \code{\link{NormMixClusK}}.
#'
#' @return
#' An S4 object of class \code{coseqResults}, with conditional
#' probabilities of cluster membership for each gene in each model stored as a list of assay
#' data, and corresponding {log likelihood, ICL value, number of
#' clusters, and form of Gaussian model} for each model stored as metadata.
#'
#' @author Andrea Rau, Cathy Maugis-Rabusseau
#'
#' @example inst/examples/NormMixClus.R
#' @importFrom BiocParallel bplapply
#' @importFrom BiocParallel bpparam
#' @export
NormMixClus <- function(y_profiles, K, subset=NULL, parallel=TRUE, BPPARAM=bpparam(),
seed = NULL, ...){
subset.index <- subset
## Parse ellipsis function
providedArgs <- list(...)
arg.user <- list(alg.type="EM", init.runs=50, init.type="small-em", init.iter=20,
iter=1000, cutoff=0.001, GaussianModel="Gaussian_pk_Lk_Ck",
verbose=TRUE, digits=3)
arg.user[names(providedArgs)] <- providedArgs
y_profiles <- as.data.frame(y_profiles)
## In case only a subset of data are to be used for analysis
if(!is.null(subset.index)) {
y_profiles <- y_profiles[subset.index,]
}
all.results <- vector("list", length = length(K))
names(all.results) <- paste0("K=", K)
if(arg.user$verbose) {
cat("Running K =", min(K), "...\n")
}
all.results[[1]] <- suppressWarnings(NormMixClusK(y_profiles, K=min(K), alg.type=arg.user$alg.type,
init.runs=arg.user$init.runs,
init.type=arg.user$init.type,
init.iter=arg.user$init.iter,
iter=arg.user$iter,
cutoff=arg.user$cutoff,
GaussianModel=arg.user$GaussianModel,
digits=arg.user$digits, seed=seed))
index <- 2
remainingK <- K[-which.min(K)]
if (length(remainingK)) {
## In the case where parallelization is NOT used
if(!parallel) {
for (k in remainingK) {
if(arg.user$verbose) {
cat("Running K =", k, "...\n")
}
all.results[[index]] <- suppressWarnings(NormMixClusK(y_profiles=y_profiles, K=k,
alg.type=arg.user$alg.type,
init.runs=arg.user$init.runs,
init.type=arg.user$init.type,
init.iter=arg.user$init.iter,
iter=arg.user$iter,
cutoff=arg.user$cutoff,
GaussianModel=arg.user$GaussianModel,
digits=arg.user$digits,
seed=seed))
index <- index + 1
}
## In the case where parallelization IS used
} else if(parallel) {
tmp <- bplapply(remainingK, function(ii, N_y_profiles, N_alg.type, N_init.runs,
N_init.type, N_init.iter, N_iter, N_cutoff,
N_GaussianModel, N_digits, N_verbose, N_FUNC, N_seed) {
if(N_verbose == TRUE) {
cat("Running K =", ii, "...\n")
}
res <- suppressWarnings(N_FUNC(y_profiles=N_y_profiles, K=as.numeric(ii),
alg.type=N_alg.type, init.runs=N_init.runs,
init.type=N_init.type, init.iter=N_init.iter,
iter=N_iter, cutoff=N_cutoff,
GaussianModel=N_GaussianModel, digits=N_digits,
seed=N_seed))
return(res)},
N_y_profiles=y_profiles, N_alg.type=arg.user$alg.type, N_init.runs = arg.user$init.runs,
N_init.type=arg.user$init.type, N_init.iter=arg.user$init.iter, N_iter=arg.user$iter,
N_cutoff=arg.user$cutoff, N_GaussianModel=arg.user$GaussianModel, N_digits=arg.user$digits,
N_verbose=arg.user$verbose, N_FUNC=NormMixClusK, N_seed = seed,
BPPARAM=BPPARAM)
if(!sum(unlist(lapply(tmp, nrow)))) {
stop(paste("All models of form", arg.user$GaussianModel, "resulted in estimation errors.
This is likely due to singular covariance matrices when this form of Gaussian mixture is used.
Try rerunning coseq with either a spherical (pk_Lk_Bk) or diagonal (pk_Lk_I)
covariance matrix form instead:
coseq(..., GaussianModel = \"Gaussian_pk_Lk_Bk\")
coseq(..., GaussianModel = \"Gaussian_pk_Lk_I\")"))
}
Kmods <- paste0("K=", unlist(lapply(tmp, function(x) nbCluster(x))))
all.results[match(Kmods, names(all.results))] <- tmp[which(unlist(lapply(tmp, function(xx)
length(nbCluster(xx)) != 0)) == TRUE)]
}
}
if(!sum(unlist(lapply(all.results, nrow)))) {
stop(paste("All models of form", arg.user$GaussianModel, "resulted in estimation errors.
This is likely due to singular covariance matrices when this form of Gaussian mixture is used.
Try rerunning coseq with either a spherical (pk_Lk_Bk) or diagonal (pk_Lk_I)
covariance matrix form instead:
coseq(..., GaussianModel = \"Gaussian_pk_Lk_Bk\")
coseq(..., GaussianModel = \"Gaussian_pk_Lk_I\")"))
}
nbClust.all <- unlist(lapply(all.results, function(x) nbCluster(x)))
logLike.all <- unlist(lapply(all.results, function(x) likelihood(x)))
ICL.all <- unlist(lapply(all.results, function(x) ICL(x)))
ICL.choose <- names(ICL.all)[which.min(ICL.all)]
select.results <- all.results[[ICL.choose]]
pp <- lapply(all.results, function(x) {
if(is.null(x)) NULL
if(!is.null(x)) assay(x)
})
names(pp) <- names(all.results)
select.results <- SummarizedExperiment(assay(select.results),
metadata=list(nbCluster=nbClust.all,
logLike=logLike.all,
ICL=ICL.all,
nbClusterError=K[!K %in% nbClust.all],
GaussianModel=arg.user$GaussianModel))
if(sum(is.nan(assay(select.results))) == prod(dim(assay(select.results)))) {
stop(paste0("The selected model, which is of the form ", arg.user$GaussianModel, ", resulted in estimation errors.
This is likely due to singular covariance matrices when this form of Gaussian mixture is used.
Try rerunning coseq with either a spherical (pk_Lk_Bk) or diagonal (pk_Lk_I)
covariance matrix form instead:
coseq(..., GaussianModel = \"Gaussian_pk_Lk_Bk\")
coseq(..., GaussianModel = \"Gaussian_pk_Lk_I\")"))
}
RES <- coseqResults(as(select.results, "RangedSummarizedExperiment"), allResults=pp)
return(RES)
}
#' Normal mixture model estimation
#'
#' Perform co-expression and co-abudance analysis of high-throughput
#' sequencing data, with or without data transformation, using a Normal
#' mixture models for single number of clusters \emph{K}.
#' The output of \code{NormMixClusK} is an S4 object of
#' class \code{RangedSummarizedExperiment}.
#'
#' @param y_profiles y (\emph{n} x \emph{q}) matrix of observed profiles for \emph{n}
#' observations and \emph{q} variables
#' @param K Number of clusters (a single value).
#' @param init.type Type of initialization strategy to be used:
#' \dQuote{\code{small-em}} for the Small-EM strategy, \dQuote{\code{random}}, \dQuote{\code{CEM}},
#' or \dQuote{\code{SEMMax}}
#' @param GaussianModel One of the 28 forms of Gaussian models defined in Rmixmod,
#' by default equal to the \code{"Gaussian_pk_Lk_Ck"} (i.e., a general family model with free
#' proportions, free volume, free shape, and free orientation)
#' @param init.runs Number of runs to be used for the Small-EM strategy, with a default value of 50
#' @param init.iter Number of iterations to be used within each run for the
#' Small-EM strategry, with a default value of 20
#' @param alg.type Algorithm to be used for parameter estimation:
#' \dQuote{\code{EM}}, \dQuote{\code{CEM}}, \dQuote{\code{SEM}}
#' @param cutoff Cutoff to declare algorithm convergence
#' @param iter Maximum number of iterations to be run for the chosen algorithm
#' @param verbose If \code{TRUE}, verbose output is created
#' @param digits Integer indicating the number of decimal places to be used for the
#' \code{probaPost} output
#' @param seed If desired, an integer defining the seed of the random number generator. If
#' \code{NULL}, a random seed is used.
#'
#' @return
#' An S4 object of class \code{RangedSummarizedExperiment}, with conditional
#' probabilities of cluster membership for each gene stored as assay data, and
#' {log likelihood, ICL value, number of
#' clusters, and form of Gaussian model} stored as metadata.
#'
#' @author Cathy Maugis-Rabusseau, Andrea Rau
#'
#' @importFrom Rmixmod mixmodGaussianModel
#' @importFrom Rmixmod mixmodStrategy
#' @importFrom Rmixmod mixmodCluster
#' @example inst/examples/NormMixClus.R
#'
#' @export
NormMixClusK <- function(y_profiles, K, alg.type="EM", init.runs=50,
init.type="small-em", GaussianModel="Gaussian_pk_Lk_Ck",
init.iter=20, iter=1000, cutoff=0.001, verbose=TRUE, digits=3,
seed = NULL) {
if(!is.data.frame(y_profiles)) y_profiles <- as.data.frame(y_profiles)
models <- mixmodGaussianModel(listModels=c(GaussianModel))
# strategy
if(init.type == "small-em") init.type <- "smallEM"
strats<-mixmodStrategy(algo=alg.type, nbTry = 1,
initMethod = init.type, nbTryInInit = init.runs,
nbIterationInInit = init.iter, nbIterationInAlgo = iter,
epsilonInInit = cutoff, epsilonInAlgo = cutoff)
if(is.null(seed)) {
xem <- mixmodCluster(y_profiles, nbCluster=K, models=models, strategy=strats,
criterion="ICL")
} else {
xem <- mixmodCluster(y_profiles, nbCluster=K, models=models, strategy=strats,
criterion="ICL", seed = seed)
}
pp <- round(proba(xem), digits=digits )
if(ncol(pp)) {
colnames(pp) <- paste0("Cluster_", seq_len(ncol(pp)))
if(!is.null(rownames(y_profiles))) rownames(pp) <- rownames(y_profiles)
if(is.null(rownames(y_profiles))) rownames(pp) <- seq_len(nrow(pp))
}
res <- SummarizedExperiment(assays=pp,
metadata=list(logLike=likelihood(xem),
ICL=ICL(xem),
nbCluster=nbCluster(xem),
GaussianModel=GaussianModel))
return(as(res, "RangedSummarizedExperiment"))
}
#' Calculate the mean and covariance for a Normal mixture model
#'
#' Calculates the mean and covariance parameters for a normal mixture model
#' of the form pK_Lk_Ck
#'
#' @param coseqResults Object of class \code{coseqResults} or \code{RangedSummarizedExperiment}
#' (as output from the \code{NormMixClus} or \code{NormMixClusK} functions)
#' @param y_profiles y (\emph{n} x \emph{q}) matrix of observed profiles for \emph{n}
#' observations and \emph{q} variables, required for \code{x} of class \code{RangedSummarizedExperiment}
#' @param K The model used for parameter estimation for objects \code{x} of
#' class \code{coseq} or \code{NormMixClus}. When \code{NULL}, the model selected
#' by the ICL criterion is used; otherwise, \code{K} should designate the number
#' of clusters in the desired model
#' @param digits Integer indicating the number of decimal places to be used for output
#' @param plot If \code{true}, produce heatmaps to visualize the estimated per-cluster
#' correlation matrices
#' @param ... Additional optional parameters to pass to \code{corrplot}, if desired
#'
#' @return
#' \item{pi }{ Vector of dimension \emph{K} with the estimated cluster proportions from
#' the Gaussian mixture model, where \emph{K} is the number of clusters}
#' \item{mu }{ Matrix of dimension \emph{K} x \emph{d} containing the estimated mean
#' vector from the Gaussian mixture model, where \emph{d} is the
#' number of samples in the data \code{y_profiles} and \emph{K} is the number of clusters}
#' \item{Sigma }{ Array of dimension \emph{d} x \emph{d} x \emph{K} containing the
#' estimated covariance matrices from the Gaussian mixture model, where \emph{d} is the
#' number of samples in the data \code{y_profiles} and \emph{K} is the number of clusters}
#' \item{rho }{ Array of dimension \emph{d} x \emph{d} x \emph{K} containing the
#' estimated correlation matrices from the Gaussian mixture model, where \emph{d} is the
#' number of samples in the data \code{y_profiles} and \emph{K} is the number of clusters}
#'
#' @author Andrea Rau, Cathy Maugis-Rabusseau
#'
#' @example inst/examples/NormMixClus.R
#' @export
#' @importFrom stats cov2cor
#' @importFrom grDevices n2mfrow
#' @importFrom graphics mtext par
NormMixParam <- function(coseqResults, y_profiles=NULL, K=NULL, digits=3, plot=FALSE, ...) {
x <- coseqResults
## TODO: object could be data.frame or matrix
if (!is(x, "coseqResults") & !is(x, "RangedSummarizedExperiment")) {
stop(paste0(sQuote("coseqResults")), " must be of class ",
paste0(dQuote("coseqResults")), " or ", paste0(dQuote("RangedSummarizedExperiment")))
}
if(is(x, "coseqResults")) {
if(nrow(tcounts(x))>0) y_profiles <- as.matrix(as.data.frame(tcounts(x)))
if(!nrow(tcounts(x)) & is.null(y_profiles)) stop("y_profiles argument needed.")
if(is.null(K)) {
probaPost <- assay(x)
} else {
if(K == "ICL") probaPost <- assay(x)
if(K != "ICL") probaPost <- coseqFullResults(x)[[paste0("K=", K)]]
}
GaussianModel <- metadata(x)$GaussianModel
}
if(is(x, "RangedSummarizedExperiment")) {
if(is.null(y_profiles)) stop(paste0(dQuote("y_profiles")), " must
be included for class ",
paste0(dQuote("RangedSummarizedExperiment")))
probaPost <- assay(x)
GaussianModel <- metadata(x)$GaussianModel
}
if(GaussianModel != "Gaussian_pk_Lk_Ck")
stop("Recalculation of Gaussian parameters currently only supported for",
paste(dQuote("Gaussian_pk_Lk_Ck")))
pi <- apply(probaPost,2,sum) / nrow(y_profiles)
mu <- matrix(0, nrow=ncol(probaPost), ncol=ncol(y_profiles))
Sigma <- array(0, dim=c(ncol(y_profiles), ncol(y_profiles), ncol(probaPost)))
for (k in seq_len(ncol(probaPost))) {
mu[k,]<-apply(probaPost[,k] * y_profiles,2,sum) / sum(probaPost[,k])
Sigma[,,k] <- (t(y_profiles) - mu[k,]) %*% ( t(t(y_profiles) - mu[k,]) * probaPost[,k])
Sigma[,,k] <- Sigma[,,k] / sum(probaPost[,k])
}
rho <- lapply(seq_len(ncol(probaPost)), function(xx) cov2cor(Sigma[,,xx]))
rho <- array(unlist(rho), dim = c(dim(rho[[1]]), length(rho)))
dimnames(rho) <- list(colnames(y_profiles), colnames(y_profiles),
paste("Cluster", seq_len(ncol(probaPost))))
dimnames(Sigma) <- list(colnames(y_profiles), colnames(y_profiles),
paste("Cluster", seq_len(ncol(probaPost))))
colnames(mu) <- colnames(y_profiles)
rownames(mu) <- paste("Cluster", seq_len(ncol(probaPost)))
if(plot) {
par(mfrow=n2mfrow(ncol(probaPost)))
for(kk in seq_len(ncol(probaPost))) {
corrplot(rho[,,kk], method="ellipse", type="upper",
tl.pos="d", ...)
mtext(paste("K =", kk), side=2, las=1, line=-1)
}
}
param <- list(pi=round(pi, digits), mu=round(mu, digits), Sigma=round(Sigma, digits),
rho=round(rho, digits))
return(param)
}
