#' Wrapper for Running Essential SCONE Modules
#'
#' @param expr matrix. The expression data matrix (genes in rows, cells in
#' columns).
#' @param qc data frame. The quality control (QC) matrix (cells in rows,
#' metrics in columns) to be used for filtering, normalization,
#' and evaluation.
#' @param bio factor. The biological condition to be modeled in the Adjustment
#' Step as variation to be preserved. If adjust_bio="no", it will not be used
#' for normalization, but only for evaluation.
#' @param batch factor. The known batch variable to be included in the
#' adjustment model as variation to be removed. If adjust_batch="no", it will
#' not be used for normalization, but only for evaluation.
#' @param negcon character. The genes to be used as negative controls for
#' filtering, normalization, and evaluation. These genes should be expressed
#' uniformily across the biological phenomenon of interest. Default NULL.
#' @param verbose character. Verbosity level: higher level is more verbose.
#' Default "0".
#' @param out_dir character. Output directory. Default getwd().
#' @param seed numeric. Random seed. Default 112233.
#'
#' @param filt_cells logical. Should cells be filtered? Set to FALSE if low
#' quality cells have already been excluded. If cells are not filtered, then
#' initial gene filtering (the one that is done prior to cell filtering) is
#' disabled as it becomes redundant with the gene filtering that is done
#' after cell filtering. Default TRUE.
#' @param filt_genes logical. Should genes be filtered post-sample filtering?
#' Default TRUE.
#' @param always_keep_genes logical. A character vector of gene names that
#' should never be excluded (e.g., marker genes). Default NULL.
#' @param fnr_maxiter numeric. Maximum number of iterations in EM estimation of
#' expression posteriors. If 0, then FNR estimation is skipped entirely, and
#' as a consequence no imputation will be performed, disregarding the value
#' of the "norm_impute" argument. Default 1000.
#'
#' @param norm_impute character. Should imputation be included in the
#' comparison? If 'force', only imputed normalizations will be run. Default
#' "yes."
#' @param norm_scaling character. Scaling options to be included in the Scaling
#' Step. Default c("none", "sum", "deseq", "tmm", "uq", "fq", "detect"). See
#' details.
#' @param norm_rezero logical. Restore prior zeroes and negative values to zero
#' following normalization. Default FALSE.
#' @param norm_k_max numeric. Max number (norm_k_max) of factors of unwanted
#' variation modeled in the Adjustment Step. Default NULL.
#' @param norm_qc_expl numeric. In automatic selection of norm_k_max, what
#' fraction of variation must be explained by the first norm_k_max PCs of qc?
#' Default 0.5. Ignored if norm_k_max is not NULL.
#' @param norm_adjust_bio character. If 'no' it will not be included in the
#' model; if 'yes', both models with and without 'bio' will be run; if
#' 'force', only models with 'bio' will be run. Default "yes."
#' @param norm_adjust_batch character. If 'no' it will not be modeled in the
#' Adjustment Step; if 'yes', both models with and without 'batch' will be
#' run; if 'force', only models with 'batch' will be run. Default "yes."
#'
#' @param eval_dim numeric. The number of principal components to use for
#' evaluation. Default NULL.
#' @param eval_expr_expl numeric. In automatic selection of eval_dim, what
#' fraction of variation must be explained by the first eval_dim PCs of expr?
#' Default 0.1. Ignored if eval_dim is not NULL.
#' @param eval_poscon character. The genes to be used as positive controls for
#' evaluation. These genes should be expected to change according to the
#' biological phenomenon of interest.
#' @param eval_max_kclust numeric. The max number of clusters (> 1) to be used
#' for pam tightness evaluation. If NULL, tightness will be returned NA.
#' @param eval_stratified_pam logical. If TRUE then maximum ASW for PAM_SIL is
#' separately computed for each biological-cross-batch condition (accepting
#' NAs), and a weighted average is returned as PAM_SIL. Default TRUE.
#' @param report_num numeric. Number of top methods to report. Default 13.
#'
#' @param out_rda logical. If TRUE, sconeResults.Rda file with the object that
#' the scone function returns is saved in the out_dir (may be very large for
#' large datasets, but useful for post-processing) Default FALSE.
#' @param eval_negcon character. Alternative negative control gene list for
#' evaluation only.
#' @param ... extra params passed to the metric_sample_filter and scone when
#' they're called by easybake
#'
#' @export
#' @importFrom grDevices dev.off pdf
#' @importFrom graphics legend lines points
#' @importFrom utils head sessionInfo write.table
#' @importFrom RColorBrewer brewer.pal
#' @importFrom boot inv.logit
#' @importFrom hexbin plot hexbin
#'
#' @details "ADD DESCRIPTION"
#'
#' @return Directory structure "ADD DESCRIPTION"
#'
#' @examples
#' set.seed(101)
#' mat <- matrix(rpois(1000, lambda = 5), ncol=10)
#' colnames(mat) <- paste("X", 1:ncol(mat), sep="")
#' obj <- SconeExperiment(mat)
#' res <- scone(obj, scaling=list(none=identity, uq=UQ_FN, deseq=DESEQ_FN),
#' evaluate=TRUE, k_ruv=0, k_qc=0, eval_kclust=2,
#' bpparam = BiocParallel::SerialParam())
#' qc = as.matrix(cbind(colSums(mat),colSums(mat > 0)))
#' rownames(qc) = colnames(mat)
#' colnames(qc) = c("NREADS","RALIGN")
#' \dontrun{
#' scone_easybake(mat, qc = as.data.frame(qc), verbose = "2",
#' norm_adjust_bio= "no",
#' norm_adjust_batch= "no", norm_k_max = 0,
#' fnr_maxiter = 0, filt_cells=FALSE, filt_genes=FALSE,
#' eval_stratified_pam = FALSE,
#' out_dir="~/scone_out")
#' }
scone_easybake <- function(expr, qc,
bio=NULL, batch=NULL, negcon=NULL,
verbose=c("0","1","2"),
out_dir=getwd(), seed = 112233,
filt_cells=TRUE,
filt_genes=TRUE, always_keep_genes = NULL,
fnr_maxiter=1000,
norm_impute=c("yes", "no", "force"),
norm_scaling=c("none", "sum", "deseq",
"tmm", "uq", "fq", "detect"),
norm_rezero=FALSE, norm_k_max=NULL, norm_qc_expl=0.5,
norm_adjust_bio=c("yes", "no", "force"),
norm_adjust_batch=c("yes", "no", "force"),
eval_dim = NULL, eval_expr_expl = 0.1,
eval_poscon = NULL, eval_negcon = negcon,
eval_max_kclust = 10, eval_stratified_pam = TRUE,
report_num = 13, out_rda = FALSE, ...) {
# Require qc or default (colSums() and colSums(>0))
verbose = match.arg(verbose)
verbose = as.numeric(verbose)
printf <- function(...) cat(sprintf(...))
set.seed(seed)
if(is.null(expr)){stop("expr must be specified")}
if(is.null(qc)){stop("qc must be specified")}
if(is.null(negcon) && fnr_maxiter > 0){
stop("negcon must be specified if fnr_maxiter > 0")
}
if(!is.matrix(expr)){stop("expr must be a matrix")}
if(!is.data.frame(qc)){stop("qc must be a data frame")}
# Primary Directory
if (!file.exists(out_dir)) {
dir.create(out_dir, recursive=TRUE)
}
# Misc Directory
misc_dir = paste0(out_dir,"/misc")
if (!file.exists(misc_dir)) {
dir.create(misc_dir)
}
#duplicate stdout to log file
logFile = file(file.path(out_dir, "stdout.txt"), open = "wt")
sink(logFile, type = "output", split = TRUE)
## ------ Data Filtering Module ------
if(filt_cells) {
# Note: There's no need to do the initial gene-filtering if
# there's no cell filtering - it becomes redundant with the
# final gene filtering (and even makes the final gene filtering
# stricter than intended because we then do gene filter twice
# on the same set of cells - increasing the quantile between
# the first and second iteration)
# Initial Gene Filtering: Select "common" transcripts.
thresh_fail = quantile(expr[expr > 0], 0.2) ###Make argument###
num_fail = 10 ###Make argument###
init.gf.vec = rowSums(expr > thresh_fail) > num_fail
if(verbose > 0){printf(paste0("Data Filtering Module: Initial filter",
" - Kept only %d genes expressed in more ",
"than %.2f units in more than %d cells, ",
"excluded %d genes\n"), sum(init.gf.vec),
thresh_fail, num_fail, sum(!init.gf.vec))}
# Initial-Filtering Negative Controls
small_negcon = intersect(negcon,rownames(expr)[init.gf.vec])
if(verbose > 0){printf(paste0("Data Filtering Module: %d negative ",
"control genes to be used for sample",
" filtering\n"), length(small_negcon))}
# Metric-based Filtering and Report
if(verbose > 0){printf("Data Filtering Module: Filtering samples\n")}
pdf(paste0(misc_dir,"/filter_report.pdf"))
mfilt = metric_sample_filter(expr, plot = TRUE, hist_breaks = 100,
zcut = 4, ###Make argument###
pos_controls = rownames(expr) %in%
small_negcon,
gene_filter = init.gf.vec,
nreads = qc$NREADS,
ralign = qc$RALIGN, ...) ###Make argument###
dev.off()
mfilt = !apply(simplify2array(mfilt[!is.na(mfilt)]),1,any)
if(verbose > 0){printf(paste0("Data Filtering Module: %d samples to be ",
"used in downstream analysis, excluded %d ",
"samples\n"), sum(mfilt), sum(!mfilt))}
# Implement sample filter
expr = expr[,mfilt]
qc = qc[mfilt,]
if(!is.null(batch)) {
batch = droplevels(batch[mfilt])
}
if(!is.null(bio)) {
bio = droplevels(bio[mfilt])
}
} else {
if(verbose > 0){printf("Data Filtering Module: Skipping cell filter...\n")}
}
if(filt_genes){
thresh_fail = quantile(expr[expr > 0], 0.2) ###Make argument###
num_fail = 10 ###Make argument###
final.gf.vec = rowSums(expr > thresh_fail) > num_fail
if(!is.null(always_keep_genes)) {
final.gf.vec[always_keep_genes %in% rownames(expr)] = TRUE
}
if(verbose > 0){printf(paste0("Data Filtering Module: Kept only %d genes ",
"expressed in more than %.2f units in more",
" than %d cells (or ones that were forced ",
"to be kept by the user's argument), ",
"excluded %d genes\n"), sum(final.gf.vec),
thresh_fail, num_fail, sum(!final.gf.vec))}
# Implement gene filter
expr = expr[final.gf.vec,]
negcon = negcon[negcon %in% rownames(expr)]
eval_negcon = eval_negcon[eval_negcon %in% rownames(expr)]
if(verbose > 0){printf(paste0("Data Filtering Module: Kept only %d ",
"negative control genes\n"), length(negcon))}
if(!is.null(eval_poscon)){
eval_poscon = eval_poscon[eval_poscon %in% rownames(expr)]
if(verbose > 0){printf(paste0("Data Filtering Module: ",
"Kept only %d positive control genes\n"),
length(eval_poscon))}
}
} else {
if(verbose > 0){printf("Data Filtering Module: Skipping gene filter...\n")}
}
if(verbose > 0){printf("Data Filtering Module: COMPLETE\n\n")}
## ------ False-Negative Rate Inference Module ------
cc <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"),brewer.pal(9,"Set3"))
if(fnr_maxiter > 0) {
## fnr_maxiter > 0 -- user wants to do fnr estimation
if(verbose > 0){printf(paste0("False-Negative Rate Inference Module:",
" Estimating posterior ",
"expression probability...\n"))}
fnr_out = estimate_ziber(x = expr,
bulk_model = TRUE, ###Make argument###
pos_controls = rownames(expr) %in%
negcon,
verbose = (verbose > 1),
maxiter = fnr_maxiter)
if(fnr_out$convergence == 1){printf(paste0("False-Negative Rate ",
"Inference Module: WARNING - ",
"EM did not converge. Try ",
"increasing fnr_max_iter\n"))}
if(verbose > 0){printf(paste0("False-Negative Rate ",
"Inference Module: Producing report...\n"))}
# FNR Report
logistic_coef = simple_FNR_params(expr = expr,
pos_controls = rownames(expr) %in%
negcon)
# "Non-DE" are expected to be expressed in all cells
pdf(paste0(misc_dir,"/fnr_report.pdf"),width = 12,height = 6)
par(mfrow=c(1,2))
plot(NULL, main = "False Negative Rate Curves",ylim = c(0,1),
xlim = c(0,15), ylab = "Failure Probability",
xlab = "Median log-Expression")
x = (0:150)/10
AUC = NULL
for(si in 1:ncol(expr)){
y = inv.logit(logistic_coef[si,1] + logistic_coef[si,2] * x)
AUC[si] = sum(y)/10
if(!is.null(batch)){
colchoice = cc[batch][si]
}else{
colchoice = "black"
}
lines(x, y, type = 'l', lwd = 2, col = colchoice)
}
# Barplot of FNR AUC
if(!is.null(batch)){
oAUC = order(AUC)
sAUC = AUC[oAUC]
sbatch = batch[oAUC]
barplot(sAUC[order(sbatch)], col=cc[sbatch][order(sbatch)],
border=cc[sbatch][order(sbatch)],
main="FNR AUC, colored by batch",
ylim = c(0,2*max(AUC)))
if(fnr_out$convergence != 1) {
# this line throws an exception
# if not converged: "'legend' is of length 0"
legend("topleft", legend=levels(sbatch), fill=cc, cex=0.4)
}
}else{
barplot(sort(AUC), col="black", border="black",
main="FNR AUC",ylim = c(0,2*max(AUC)))
}
hexbin::plot(hexbin(rowMeans(expr > 0)[!is.na(fnr_out$Beta)],
inv.logit(fnr_out$Beta[!is.na(fnr_out$Beta)])),
xlab = "Detection Rate", ylab = "Expression Rate")
hexbin::plot(hexbin(fnr_out$Alpha[1,][!is.na(fnr_out$Beta)],
inv.logit(fnr_out$Beta[!is.na(fnr_out$Beta)])),
ylab = "Expression Rate", xlab = "Median log-Expression")
if(verbose > 0){printf(paste0("False-Negative Rate ",
"Inference Module: COMPLETE\n\n"))}
dev.off()
} else {
## fnr_maxiter <= 0 (should be 0, check for <= 0 for robustness)
## -- user doesn't want to do fnr estimation
if(verbose > 0) {
printf("Skipping False-Negative Rate Inference Module...\n")
printf("(and consequently no imputation will be performed too)\n\n")
norm_impute = "no"
}
}
## ------ SCONE Prime Module ------
qcmat = as.matrix(qc)
#Allow non-serial (can register serial from outside for debug)
#BiocParallel::register(BiocParallel::SerialParam()) ###Make argument###
# Imputation arguments
norm_impute = match.arg(norm_impute)
imputation = switch(norm_impute,
yes=list(none=impute_null,expect=impute_expectation),
no=list(none=impute_null),
force=list(expect=impute_expectation))
if((norm_impute %in% c("yes", "force")) && fnr_maxiter > 0) {
impute_args = list(p_nodrop = fnr_out$p_nodrop,
mu = exp(fnr_out$Alpha[1,]))
} else {
impute_args = NULL
}
# Scaling arguments
norm_scaling = match.arg(norm_scaling,several.ok = TRUE)
SUM_FN = function (ei)
{
sums = colSums(ei)
eo = t(t(ei)*mean(sums)/sums)
return(eo)
}
EFF_FN = function (ei)
{
sums = colSums(ei > 0)
eo = t(t(ei)*sums/mean(sums))
return(eo)
}
scaling=list(none=identity,
sum = SUM_FN,
deseq=DESEQ_FN,
tmm = TMM_FN,
uq = UQ_FN,
fq = FQT_FN,
detect = EFF_FN)[norm_scaling]
# Adjustment arguments
if(is.null(norm_k_max)){
qc_sd = prcomp(qcmat,center = TRUE,scale = TRUE)$sd
norm_k_max = which(cumsum(qc_sd^2)/sum(qc_sd^2) > norm_qc_expl)[1]
}
if(verbose > 0){printf(paste0("Normalization Module: Adjusting ",
"for up to %d factors",
" of unwanted variation\n"),
norm_k_max)}
# Generate Params (RUN = FALSE)
if(verbose > 0){printf("Normalization Module: Selecting params:\n")}
args_list = list(object = expr,
qc=qcmat)
# Creating a SconeExperiment Object
if(!is.null(negcon)){
args_list = c( args_list,
list( negcon_ruv = rownames(expr) %in% negcon ))
}
if(!is.null( eval_negcon)){
args_list = c( args_list,
list( negcon_eval = rownames(expr) %in% eval_negcon ))
}
if(!is.null(eval_poscon)){
args_list = c( args_list,
list( poscon = rownames(expr)%in% eval_poscon ))
}
if(!is.null(bio)){
args_list = c( args_list, list( bio = bio ))
}
if(!is.null(batch)){
args_list = c( args_list, list( batch = batch ))
}
my_scone <- do.call(SconeExperiment,args_list)
my_scone <- scone(my_scone,
imputation = imputation, impute_args = impute_args,
scaling=scaling,
k_qc=norm_k_max, k_ruv = norm_k_max,
adjust_bio = match.arg(norm_adjust_bio),
adjust_batch = match.arg(norm_adjust_batch),
eval_kclust = NULL,
run=FALSE, ...)
is_screened = ((get_params(my_scone)$imputation_method == "expect") &
(get_params(my_scone)$scaling_method %in% c("detect"))) |
((get_params(my_scone)$adjust_biology == "bio")
& (get_params(my_scone)$adjust_batch != "batch"))
if(norm_rezero){
zerop = "postadjust"
}else{
zerop = "none"
}
my_scone = select_methods(my_scone,
rownames(get_params(my_scone))[!is_screened ])
if(verbose > 0){print(get_params(my_scone))}
# Evaluation arguments
if(is.null(eval_dim)){
expr_sd = prcomp(t(expr),center = TRUE,scale = TRUE)$sd
eval_dim = which(cumsum(expr_sd^2)/sum(expr_sd^2) > eval_expr_expl)[1]
}
if(verbose > 0){printf(paste0("Normalization Module: Evaluation ",
"based on first %d PCs of expression\n"),
eval_dim)}
if(is.null(eval_max_kclust)){
eval_kclust = NULL
if(verbose > 0){printf("Normalization Module: No PAM-based evaluation\n")}
}else{
if(eval_stratified_pam){
if(!is.null(bio) && !is.null(batch)) {
temp_lim = min(table(bio,batch))
} else if(is.null(batch) && is.null(bio)) {
temp_lim = ncol(expr)
} else if (!is.null(bio)) { #one of them is not null and the other is
temp_lim = min(table(bio))
} else {
temp_lim = min(table(batch))
}
uplim = min(temp_lim-1,eval_max_kclust)
}else{
uplim = min(ncol(expr)-1,eval_max_kclust)
}
if(uplim >=2){
eval_kclust = 2:uplim
if(verbose > 0){printf(paste0("Normalization Module: ",
" Maximum clustering k for",
" PAM set to %d\n"),
uplim)}
}else{
eval_kclust = NULL
if(verbose > 0){printf(paste0("Normalization Module:",
" No PAM-based evaluation\n"))}
}
}
# Generate Scores and Ranking
if(verbose > 0){printf("Normalization Module: Scoring Normalizations...\n")}
tic = proc.time()
my_scone <- scone(my_scone,
imputation = imputation, impute_args = impute_args,
scaling=scaling,
k_qc=norm_k_max, k_ruv = norm_k_max,
adjust_bio = match.arg(norm_adjust_bio),
adjust_batch = match.arg(norm_adjust_batch),
run=TRUE, verbose = (verbose > 1),
stratified_pam = eval_stratified_pam,
eval_kclust = eval_kclust,
zero = zerop, ...)
toc = proc.time()
if(verbose > 0) {
printf("Normalization Module: Scoring Normalizations Done!...\n")
printf("time elapsed (seconds):\n")
print(toc-tic)
}
# Generate SCONE Report
if(verbose > 0){printf("Normalization Module: Producing main report...\n")}
pdf(paste0(out_dir,"/scone_report.pdf"),width = 6,height = 6)
pc_obj = prcomp(apply(t(get_scores(my_scone)),1,rank),
center = TRUE,scale = FALSE)
bp_obj = biplot_color(pc_obj,y = -get_score_ranks(my_scone),expand = .6)
points(t(bp_obj[1,]), pch = 1, col = "red", cex = 1)
points(t(bp_obj[1,]), pch = 1, col = "red", cex = 1.5)
points(t(bp_obj[rownames(bp_obj) == "none,none,no_uv,no_bio,no_batch",]),
pch = 1, col = "blue", cex = 1)
points(t(bp_obj[rownames(bp_obj) == "none,none,no_uv,no_bio,no_batch",]),
pch = 1, col = "blue", cex = 1.5)
arrows(bp_obj[rownames(bp_obj) == "none,none,no_uv,no_bio,no_batch",][1],
bp_obj[rownames(bp_obj) == "none,none,no_uv,no_bio,no_batch",][2],
bp_obj[1,][1],
bp_obj[1,][2],
lty = 2, lwd = 2)
dev.off()
# Recomputing top methods
params_select = head(get_params(my_scone),n = report_num)
if(verbose > 0){printf(paste0("Normalization Module: ",
"Re-computing top normalizations...\n"))}
tic = proc.time()
scone_res = list()
# List of normalized matrices
scone_res$normalized_data = lapply(as.list(rownames(params_select)),
FUN = function(z){
get_normalized(my_scone,
method = z,log=TRUE)
})
names(scone_res$normalized_data) = rownames(params_select)
# Parameter matrix
scone_res$params = get_params(my_scone)[rownames(params_select),]
# Merged score matrix
scone_res$scores = cbind(get_scores(my_scone),
get_score_ranks(my_scone))[rownames(params_select),]
colnames(scone_res$scores)[ncol(scone_res$scores)] = "mean_score_rank"
toc = proc.time()
if(verbose > 0) {
printf("Normalization Module: Recomputing top methods done!...\n")
printf("time elapsed (seconds):\n")
print(toc-tic)
}
if(out_rda) {
if(verbose > 0){printf(paste0("Normalization Module: ",
"Writing sconeResults.Rda file ",
"of scone's returned object...\n"))}
save(file=file.path(out_dir, "sconeResults.Rda"), scone_res)
}
if(verbose > 0){
printf("Normalization Module: (Trumpets): the selected method is... %s!\n",
rownames( scone_res$scores)[1])
printf("Normalization Module: Writing top normalization to file...\n")
}
normle = scone_res$normalized_data[[1]]
write.table(x = normle,row.names = TRUE,col.names = TRUE,quote = FALSE,
sep = "\t",eol = "\n",file = paste0(out_dir,"/best_scone.txt"))
if(verbose > 0){printf(paste0("Normalization Module: ",
"Producing normalization-specific ",
"reports...\n"))}
write.table(x = scone_res$scores, row.names = TRUE, col.names = TRUE,
quote = FALSE, sep = "\t", eol = "\n",
file = file.path(out_dir,"/normalization_scores.txt"))
for(count in seq_along(rownames( scone_res$scores))) {
nom = rownames( scone_res$scores)[count]
if(verbose > 0){
printf(paste0("Normalization Module:\t(",count,")\t",nom,"\n"))
}
nom_dir = paste0(misc_dir,"/N",count,"_",gsub(",","_",nom))
if (!file.exists(nom_dir)) {
dir.create(nom_dir)
}
normle = scone_res$normalized_data[[nom]]
write.table(x = normle,row.names = TRUE,col.names = TRUE,
quote = FALSE,sep = "\t",eol = "\n",
file = paste0(nom_dir,"/norm.txt"))
if(is.null(bio) & is.null(batch)){
pdf(paste0(nom_dir,"/norm_report.pdf"),width = 6,height = 6)
plot(prcomp(t(normle ))$x, col = "black", pch =16)
dev.off()
}else if(is.null(bio) & !is.null(batch)){
pdf(paste0(nom_dir,"/norm_report.pdf"),width = 6,height = 6)
plot(prcomp(t(normle ))$x, col = cc[batch], pch =16)
dev.off()
}else if(!is.null(bio) & is.null(batch)){
pdf(paste0(nom_dir,"/norm_report.pdf"),width = 6,height = 6)
plot(prcomp(t(normle ))$x, col = cc[bio], pch =16)
dev.off()
}else{
pdf(paste0(nom_dir,"/norm_report.pdf"),width = 12,height = 6)
par(mfrow=c(1,2))
plot(prcomp(t(normle ))$x, col = cc[batch], pch =16)
plot(prcomp(t(normle ))$x, col = cc[bio], pch =16)
dev.off()
}
}
if(verbose > 0){printf("Normalization Module: COMPLETE\n")}
## ------ Bye now ------
if(verbose > 0){
printf("\n=====sessionInfo=====\n\n")
sessionInfo()
}
#stop duplicating stdout to log file
sink(type = "output")
#sink(type = "message") -- no, cannot split the message connection
}
