@@ -1,5 +1,5 @@

 Package: scone

-Version: 1.1.2

+Version: 1.1.3

 Title: Single Cell Overview of Normalized Expression data

 Description: SCONE is an R package for comparing and ranking the performance of

 	different normalization schemes for single-cell RNA-seq and other 

@@ -60,6 +60,7 @@ importFrom(grDevices,pdf)

 importFrom(graphics,abline)

 importFrom(graphics,arrows)

 importFrom(graphics,barplot)

+importFrom(graphics,boxplot)

 importFrom(graphics,hist)

 importFrom(graphics,legend)

 importFrom(graphics,lines)

@@ -91,6 +92,7 @@ importFrom(stats,dist)

 importFrom(stats,dnbinom)

 importFrom(stats,fitted.values)

 importFrom(stats,glm)

+importFrom(stats,lm)

 importFrom(stats,mad)

 importFrom(stats,median)

 importFrom(stats,model.matrix)

@@ -101,7 +103,10 @@ importFrom(stats,prcomp)

 importFrom(stats,quantile)

 importFrom(stats,quasibinomial)

 importFrom(stats,sd)

+importFrom(stats,setNames)

+importFrom(stats,var)

 importFrom(utils,capture.output)

 importFrom(utils,head)

 importFrom(utils,sessionInfo)

+importFrom(utils,write.csv)

 importFrom(utils,write.table)

@@ -33,7 +33,8 @@ TMM_FN = function(ei) {

   return(eo)

 }

-#' Relative log-expression (RLE; DESeq) scaling normalization wrapper function

+#' Relative log-expression (RLE; DESeq) scaling normalization wrapper

+#' function

 #' @importFrom edgeR calcNormFactors

 #' @details SCONE scaling wrapper for \code{\link[edgeR]{calcNormFactors}}).

 #' @export

@@ -24,6 +24,9 @@

 #'

 #' @importFrom RColorBrewer brewer.pal

 #' @importFrom rARPACK svds

+#' @importFrom graphics boxplot

+#' @importFrom utils write.csv

+#' @importFrom stats setNames var

 #' @export

 #'

 #' @return An object that represents the SCONE report app.

@@ -66,15 +69,15 @@ sconeReport = function(x, methods,

     stop("reshape2 package needed for sconeReport()")

   }

-  if (!require("plotly", quietly = TRUE)) {

+  if (!requireNamespace("plotly", quietly = TRUE)) {

     stop("plotly package needed for sconeReport()")

   }

-  if (!require("visNetwork", quietly = TRUE)) {

+  if (!requireNamespace("visNetwork", quietly = TRUE)) {

     stop("visNetwork package needed for sconeReport()")

   }

-  if (!require("ggplot2", quietly = TRUE)) {

+  if (!requireNamespace("ggplot2", quietly = TRUE)) {

     stop("ggplot2 package needed for sconeReport()")

   }

@@ -223,18 +226,18 @@ sconeReport = function(x, methods,

                             ),

                             shiny::tabPanel("PCA",

                                      shiny::br(),

-                                     shiny::p(paste0("This panel shows principal ",

+                                shiny::p(paste0("This panel shows principal ",

                                               "component analysis ",

                                               "(PCA) results",

                                               " on different subsets ",

                                               "of genes.")),

                                      shiny::br(),

-                                     shiny::h6("Variance Explained (All Genes)"),

+                                  shiny::h6("Variance Explained (All Genes)"),

                                      shiny::p(paste0("Use this plot",

                                               " to decide on the ",

                                               "dimensionality of the reduced",

                                               " spaced used for evaluation.")),

-                                     shiny::plotOutput("plot_scree",width = "650px",

+                                shiny::plotOutput("plot_scree",width = "650px",

                                                 height = "400px"),

                                      shiny::br(),

                                      shiny::h6("2D (All Genes)"),

@@ -264,7 +267,7 @@ sconeReport = function(x, methods,

                             ),

                             shiny::tabPanel("QC",

                                      shiny::br(),

-                                     shiny::p(paste0("This panel shows the absolute",

+                              shiny::p(paste0("This panel shows the absolute",

                                               " correlations between",

                                               " Principal",

                                               " Components (PCs) of",

@@ -326,7 +329,7 @@ sconeReport = function(x, methods,

                             ),

                             shiny::tabPanel("Control Genes",

                                      shiny::br(),

-                                     shiny::p(paste0("Heatmap of control genes, ",

+                                  shiny::p(paste0("Heatmap of control genes, ",

                                               "colored by all",

                                               " three categories.")),

                                      shiny::br(),

@@ -625,7 +628,7 @@ sconeReport = function(x, methods,

              fill = cc[sort(unique(factor(strat_col())))])

     })

-    output$plot3d_base <- shiny::renderPlotly({

+    output$plot3d_base <- plotly::renderPlotly({

       PC1 <- PC2 <- PC3 <- NULL

       df <- setNames(data.frame(pc_obj_base()$x[,1:3]),

                      c("PC1", "PC2", "PC3"))

@@ -665,7 +668,7 @@ sconeReport = function(x, methods,

       abs(cor(qc,pc_obj_base()$x))

     })

-    output$qccorPlot <- shiny::renderPlotly({

+    output$qccorPlot <- plotly::renderPlotly({

       metric <- value <- PC <- NULL

       df = reshape2::melt(cor_qc()[,1:input$dim])

@@ -696,7 +699,7 @@ sconeReport = function(x, methods,

              fill = cc[sort(unique(factor(strat_col())))])

     })

-    output$plot3d_qc <- shiny::renderPlotly({

+    output$plot3d_qc <- plotly::renderPlotly({

       if(ncol(pc_obj_qc()$x) >= 3){

         PC1 <- PC2 <- PC3 <- NULL

         df <- setNames(data.frame(pc_obj_qc()$x[,1:3]),

@@ -1,120 +1,146 @@

 #' SCONE Evaluation: Evaluate an Expression Matrix

-#' 

-#' This function evaluates a (normalized) expression matrix using SCONE 

-#' criteria, producing 8 metrics based on i) Clustering, ii) Correlations and 

+#'

+#' This function evaluates a (normalized) expression matrix using SCONE

+#' criteria, producing 8 metrics based on i) Clustering, ii) Correlations and

 #' iii) Relative Expression.

-#' 

+#'

 #' @details Users may specify their own eval_proj function that will be used to

 #'   compute Clustering and Correlation metrics. This eval_proj() function must

-#'   have 2 input arguments: \itemize{ \item{e}{ matrix. log-transformed (+ 

-#'   pseudocount) expression data (genes in rows, cells in columns).} 

+#'   have 2 input arguments: \itemize{ \item{e}{ matrix. log-transformed (+

+#'   pseudocount) expression data (genes in rows, cells in columns).}

 #'   \item{eval_proj_args}{ list. additional function arguments, e.g. prior

 #'   data weights.}} and it must output a matrix representation of the original

 #'   data (cells in rows, factors in columns). The value of eval_proj_args is

 #'   passed to the user-defined function from the eval_proj_args argument of

 #'   the main score_matrix() function call.

-#'   

-#' @param expr matrix. The expression data matrix (genes in rows, cells in 

+#'

+#' @param expr matrix. The expression data matrix (genes in rows, cells in

 #'   columns).

 #' @param eval_pcs numeric. The number of principal components to use for

 #'   evaluation (Default 3). Ignored if !is.null(eval_proj).

 #' @param eval_proj function. Projection function for evaluation (see Details).

 #'   If NULL, PCA is used for projection

-#' @param eval_proj_args list. List of arguments passed to projection function 

+#' @param eval_proj_args list. List of arguments passed to projection function

 #'   as eval_proj_args (see Details).

-#' @param eval_kclust numeric. The number of clusters (> 1) to be used for pam 

-#'   tightness (PAM_SIL) evaluation. If an array of integers, largest average 

+#' @param eval_kclust numeric. The number of clusters (> 1) to be used for pam

+#'   tightness (PAM_SIL) evaluation. If an array of integers, largest average

 #'   silhouette width (tightness) will be reported in PAM_SIL. If NULL, PAM_SIL

 #'   will be returned NA.

 #' @param bio factor. A known biological condition (variation to be preserved),

 #'   NA is allowed. If NULL, condition ASW, BIO_SIL, will be returned NA.

 #' @param batch factor. A known batch variable (variation to be removed), NA is

 #'   allowed. If NULL, batch ASW, BATCH_SIL, will be returned NA.

-#' @param qc_factors Factors of unwanted variation derived from quality 

+#' @param qc_factors Factors of unwanted variation derived from quality

 #'   metrics. If NULL, qc correlations, EXP_QC_COR, will be returned NA.

-#' @param uv_factors Factors of unwanted variation derived from negative 

+#' @param uv_factors Factors of unwanted variation derived from negative

 #'   control genes (evaluation set). If NULL, uv correlations, EXP_UV_COR,

 #'   will be returned NA.

 #' @param wv_factors Factors of wanted variation derived from positive

 #'   control genes (evaluation set). If NULL, wv correlations, EXP_WV_COR,

 #'   will be returned NA.

-#' @param is_log logical. If TRUE the expr matrix is already logged and log 

+#' @param is_log logical. If TRUE the expr matrix is already logged and log

 #'   transformation will not be carried out prior to projection. Default FALSE.

-#' @param stratified_pam logical. If TRUE then maximum ASW is separately 

-#'   computed for each biological-cross-batch stratum (accepts NAs), and a 

+#' @param stratified_pam logical. If TRUE then maximum ASW is separately

+#'   computed for each biological-cross-batch stratum (accepts NAs), and a

 #'   weighted average silhouette width is returned as PAM_SIL. Default FALSE.

-#'   

+#' @param stratified_cor logical. If TRUE then cor metrics are separately

+#'   computed for each biological-cross-batch stratum (accepts NAs), and

+#'   weighted averages are returned for EXP_QC_COR, EXP_UV_COR, & EXP_WV_COR.

+#'   Default FALSE.

+#' @param stratified_rle logical. If TRUE then rle metrics are separately

+#'   computed for each biological-cross-batch stratum (accepts NAs), and

+#'   weighted averages are returned for RLE_MED & RLE_IQR. Default FALSE.

+#'

 #' @importFrom class knn

 #' @importFrom fpc pamk

 #' @importFrom cluster silhouette

 #' @importFrom rARPACK svds

 #' @importFrom matrixStats rowMedians colMedians colIQRs

-#'   

+#' @importFrom stats lm

+#'

 #' @export

-#' 

+#'

 #' @return A list with the following metrics: \itemize{ \item{BIO_SIL}{ Average

-#'   silhouette width by biological condition.} \item{BATCH_SIL}{ Average 

-#'   silhouette width by batch condition.} \item{PAM_SIL}{ Maximum average 

-#'   silhouette width from PAM clustering (see stratified_pam argument).} 

-#'   \item{EXP_QC_COR}{ Maximum squared spearman correlation between expression

-#'   pcs and quality factors.} \item{EXP_UV_COR}{ Maximum squared spearman 

-#'   correlation between expression pcs and negative control gene factors.} 

-#'   \item{EXP_WV_COR}{ Maximum squared spearman correlation between expression

-#'   pcs and positive control gene factors.} \item{RLE_MED}{ The mean squared 

-#'   median Relative Log Expression (RLE).} \item{RLE_IQR}{ The variance of the

-#'   inter-quartile range (IQR) of the RLE.} }

-#'   

-#' @examples 

-#' 

+#'   silhouette width by biological condition.} \item{BATCH_SIL}{ Average

+#'   silhouette width by batch condition.} \item{PAM_SIL}{ Maximum average

+#'   silhouette width from PAM clustering (see stratified_pam argument).}

+#'   \item{EXP_QC_COR}{ Coefficient of determination between expression

+#'   pcs and quality factors (see stratified_cor argument).} \item{EXP_UV_COR}{

+#'   Coefficient of determination between expression pcs and negative

+#'   control gene factors (see stratified_cor argument).} \item{EXP_WV_COR}{

+#'   Coefficient of determination between expression pcs and positive

+#'   control gene factors (see stratified_cor argument).} \item{RLE_MED}{ The

+#'   mean squared median Relative Log Expression (RLE) (see stratified_rle

+#'   argument).} \item{RLE_IQR}{ The variance of the inter-quartile range (IQR)

+#'   of the RLE (see stratified_rle argument).} }

+#'

+#' @examples

+#'

 #' set.seed(141)

 #' bio = as.factor(rep(c(1,2),each = 2))

 #' batch = as.factor(rep(c(1,2),2))

 #' log_expr = matrix(rnorm(20),ncol = 4)

-#' 

-#' scone_metrics = score_matrix(log_expr, 

+#'

+#' scone_metrics = score_matrix(log_expr,

 #'    bio = bio, batch = batch,

 #'    eval_kclust = 2, is_log = TRUE)

-#' 

+#'

-score_matrix <- function(expr, eval_pcs = 3,

-                         eval_proj = NULL, eval_proj_args = NULL,

+score_matrix <- function(expr,

+                         eval_pcs = 3,

+                         eval_proj = NULL,

+                         eval_proj_args = NULL,

                          eval_kclust = NULL,

-                         bio = NULL, batch = NULL,

+                         bio = NULL,

+                         batch = NULL,

                          qc_factors = NULL,

                          uv_factors = NULL,

                          wv_factors = NULL,

-                         is_log=FALSE, stratified_pam = FALSE){

-  

-  if(any(is.na(expr) | is.infinite(expr) | is.nan(expr))){

+                         is_log = FALSE,

+                         stratified_pam = FALSE,

+                         stratified_cor = FALSE,

+                         stratified_rle = FALSE) {

+  if (any(is.na(expr) | is.infinite(expr) | is.nan(expr))) {

     stop("NA/Inf/NaN Expression Values.")

   }

-  if(!is_log) {

+  if (!is_log) {

     expr <- log1p(expr)

   }

-  # The svd we do below on expr throws an exception if 

-  # expr created by one of the normalizations has a 

+  # The svd we do below on expr throws an exception if

+  # expr created by one of the normalizations has a

   # constant feature (=gene, i.e. row)

-  constantFeatures = apply(expr, 1, function(x) max(x)-min(x)) < 1e-3

-  if(any(constantFeatures)) {

-    warning(sprintf(paste0("scone_eval: expression matrix ",

-                           "contained %d constant features (rows) ",

-                           "---> excluding them"), sum(constantFeatures)))

-    expr = expr[!constantFeatures, ]

+  constantFeatures = apply(expr, 1, function(x)

+    max(x) - min(x)) < 1e-3

+  if (any(constantFeatures)) {

+    warning(sprintf(

+      paste0(

+        "scone_eval: expression matrix ",

+        "contained %d constant features (rows) ",

+        "---> excluding them"

+      ),

+      sum(constantFeatures)

+    ))

+    expr = expr[!constantFeatures,]

   }

-  if(is.null(eval_proj)){

+  if (is.null(eval_proj)) {

     proj = tryCatch({

-      svds(scale(t(expr), center = TRUE, scale = TRUE),

-           k=eval_pcs, nu=eval_pcs, nv=0)$u},

-      error = function(e) {

-        stop("scone_eval: svd failed")

-      })

+      svds(

+        scale(t(expr), center = TRUE, scale = TRUE),

+        k = eval_pcs,

+        nu = eval_pcs,

+        nv = 0

+      )$u

+    },

+    error = function(e) {

+      stop("scone_eval: svd failed")

+    })

   } else {

-    proj = eval_proj(expr,eval_proj_args = eval_proj_args)

+    proj = eval_proj(expr, eval_proj_args = eval_proj_args)

     eval_pcs = ncol(proj)

   }

@@ -123,17 +149,21 @@ score_matrix <- function(expr, eval_pcs = 3,

   # Biological Condition

-  if( !is.null(bio) ) {

-    if(!all(is.na(bio))) {

-      if(length(unique(bio)) > 1) {

+  if (!is.null(bio)) {

+    if (!all(is.na(bio))) {

+      if (length(unique(bio)) > 1) {

         BIO_SIL = summary(cluster::silhouette(as.numeric(na.omit(bio)),

                                               dd[!is.na(bio),

                                                  !is.na(bio)]))$avg.width

       } else {

         BIO_SIL = NA

-        warning(paste0("after exclusion of samples, ",

-                       "only one bio remains, BATCH_BIO",

-                       " is undefined"))

+        warning(

+          paste0(

+            "after exclusion of samples, ",

+            "only one bio remains, BATCH_BIO",

+            " is undefined"

+          )

+        )

       }

     } else {

       BIO_SIL = NA

@@ -145,17 +175,21 @@ score_matrix <- function(expr, eval_pcs = 3,

   # Batch Condition

-  if(!is.null(batch)) {

-    if(!all(is.na(batch))) {

-      if(length(unique(batch)) > 1) {

+  if (!is.null(batch)) {

+    if (!all(is.na(batch))) {

+      if (length(unique(batch)) > 1) {

         BATCH_SIL <- summary(cluster::silhouette(as.numeric(na.omit(batch)),

                                                  dd[!is.na(batch),

                                                     !is.na(batch)]))$avg.width

       } else {

         BATCH_SIL <- NA

-        warning(paste0("after exclusion of samples,",

-                       " only one batch remains, ",

-                       "BATCH_SIL is undefined"))

+        warning(

+          paste0(

+            "after exclusion of samples,",

+            " only one batch remains, ",

+            "BATCH_SIL is undefined"

+          )

+        )

       }

     } else{

       BATCH_SIL <- NA

@@ -167,94 +201,197 @@ score_matrix <- function(expr, eval_pcs = 3,

   ## ------ PAM Tightness -----

-  if ( !is.null(eval_kclust) ){

-    

+  if (!is.null(eval_kclust)) {

     # "Stratified" PAM

-    if(stratified_pam){

-      

-      biobatch = as.factor(paste(bio,batch,sep = "_"))

+    if (stratified_pam) {

+      biobatch = as.factor(paste(bio, batch, sep = "_"))

       PAM_SIL = 0

       # Max Average Sil Width per bio-cross-batch Condition

-      for (cond in levels(biobatch)){

+      for (cond in levels(biobatch)) {

         is_cond = which(biobatch == cond)

         cond_w = length(is_cond)

-        if(cond_w > max(eval_kclust)){

-          pamk_object = pamk(proj[is_cond,],krange = eval_kclust)

+        if (cond_w > max(eval_kclust)) {

+          pamk_object = pamk(proj[is_cond, ], krange = eval_kclust)

-          # Despite krange excluding nc = 1, 

+          # Despite krange excluding nc = 1,

           # if asw is negative, nc = 1 will be selected

-          if(is.null(pamk_object$pamobject$silinfo$avg.width) ){

+          if (is.null(pamk_object$pamobject$silinfo$avg.width)) {

             if (!1 %in% eval_kclust) {

-              pamk_object$pamobject$silinfo$avg.width = max(

-                pamk_object$crit[1:max(eval_kclust) %in% eval_kclust]

-              )

-            }else{

-              stop(paste0("nc = 1 was selected by Duda-Hart,",

-                          " exclude 1 from eval_kclust."))

+              pamk_object$pamobject$silinfo$avg.width = 

+                max(pamk_object$crit[1:max(eval_kclust) %in% eval_kclust])

+            } else{

+              stop(paste0(

+                "nc = 1 was selected by Duda-Hart,",

+                " exclude 1 from eval_kclust."

+              ))

             }

           }

-          PAM_SIL = PAM_SIL + cond_w*pamk_object$pamobject$silinfo$avg.width

-        }else{

-          stop(paste(paste0("Number of clusters 'k' must be ",

-                            "smaller than bio-cross-batch stratum size:"),

-                     paste(levels(biobatch),

-                           table(biobatch),

-                           sep = " = ",collapse = ", ")))

+          PAM_SIL = PAM_SIL + cond_w * pamk_object$pamobject$silinfo$avg.width

+        } else{

+          stop(paste(

+            paste0(

+              "Number of clusters 'k' must be ",

+              "smaller than bio-cross-batch stratum size:"

+            ),

+            paste(

+              levels(biobatch),

+              table(biobatch),

+              sep = " = ",

+              collapse = ", "

+            )

+          ))

         }

       }

-      PAM_SIL = PAM_SIL/length(biobatch)

+      PAM_SIL = PAM_SIL / length(biobatch)

       # Traditional PAM

-    }else{

-      pamk_object = pamk(proj,krange = eval_kclust)

+    } else{

+      pamk_object = pamk(proj, krange = eval_kclust)

       PAM_SIL = pamk_object$pamobject$silinfo$avg.width

     }

-  }else{

+  } else{

     PAM_SIL = NA

   }

   ## ------ Correlation with Factors -----

-  

-  # Max cor with quality factors.

-  if(!is.null(qc_factors)){

-    EXP_QC_COR = max(cor(proj,qc_factors,method = "spearman")^2)

-  }else{

-    EXP_QC_COR = NA

-  }

-  

-  # Max cor with UV factors.

-  if(!is.null(uv_factors)){

-    EXP_UV_COR = max(cor(proj,uv_factors,method = "spearman")^2)

-  }else{

-    EXP_UV_COR = NA

-  }

-  

-  # Max cor with WV factors.

-  if(!is.null(wv_factors)){

-    EXP_WV_COR = max(cor(proj,wv_factors,method = "spearman")^2)

-  }else{

-    EXP_WV_COR = NA

+  if (stratified_cor) {

+    biobatch = as.factor(paste(bio, batch, sep = "_"))

+    EXP_QC_COR <- EXP_UV_COR <- EXP_WV_COR <- 0

+    

+    for (cond in levels(biobatch)) {

+      is_cond = which(biobatch == cond)

+      cond_w = length(is_cond)

+      

+      # Max cor with quality factors.

+      if (!is.null(qc_factors)) {

+        EXP_QC_COR <-

+          EXP_QC_COR + cond_w * (1 - sum(unlist(apply(proj[is_cond,], 2,

+                                                      function(y) {

+            lm(y ~ qc_factors[is_cond,])$residual

+          })) ^ 2) / sum(scale(proj[is_cond,], scale = FALSE) ^ 2))

+      }

+      

+      # Max cor with UV factors.

+      if (!is.null(uv_factors)) {

+        EXP_UV_COR  <-

+          EXP_UV_COR + cond_w * (1 - sum(unlist(apply(proj[is_cond,], 2,

+                                                      function(y) {

+            lm(y ~ uv_factors[is_cond,])$residual

+          })) ^ 2) / sum(scale(proj[is_cond,], scale = FALSE) ^ 2))

+      }

+      

+      # Max cor with WV factors.

+      if (!is.null(wv_factors)) {

+        EXP_WV_COR <-

+          EXP_WV_COR + cond_w * (1 - sum(unlist(apply(proj[is_cond,], 2,

+                                                      function(y) {

+            lm(y ~ wv_factors[is_cond,])$residual

+          })) ^ 2) / sum(scale(proj[is_cond,], scale = FALSE) ^ 2))

+      }

+      

+    }

+    

+    if (!is.null(qc_factors)) {

+      EXP_QC_COR <- EXP_QC_COR / length(biobatch)

+    } else{

+      EXP_QC_COR <- NA

+    }

+    if (!is.null(uv_factors)) {

+      EXP_UV_COR <- EXP_UV_COR / length(biobatch)

+    } else{

+      EXP_UV_COR <- NA

+    }

+    if (!is.null(wv_factors)) {

+      EXP_WV_COR <- EXP_WV_COR / length(biobatch)

+    } else{

+      EXP_WV_COR <- NA

+    }

+    

+  } else{

+    # Max cor with quality factors.

+    if (!is.null(qc_factors)) {

+      EXP_QC_COR <- 1 - sum(unlist(apply(proj, 2, function(y) {

+        lm(y ~ qc_factors)$residual

+        })) ^ 2) / sum(scale(proj, scale = FALSE) ^ 2)

+    } else{

+      EXP_QC_COR = NA

+    }

+    

+    # Max cor with UV factors.

+    if (!is.null(uv_factors)) {

+      EXP_UV_COR <- 1 - sum(unlist(apply(proj, 2, function(y) {

+        lm(y ~ uv_factors)$residual

+      })) ^ 2) / sum(scale(proj, scale = FALSE) ^ 2)

+    } else{

+      EXP_UV_COR = NA

+    }

+    

+    # Max cor with WV factors.

+    if (!is.null(wv_factors)) {

+      EXP_WV_COR <- 1 - sum(unlist(apply(proj, 2, function(y) {

+        lm(y ~ wv_factors)$residual

+      })) ^ 2) / sum(scale(proj, scale = FALSE) ^ 2)

+    } else{

+      EXP_WV_COR = NA

+    }

+    

   }

   ## ----- RLE Measures

-  rle <- expr - rowMedians(expr)

-  

-  # Non-Zero Median RLE

-  RLE_MED <- mean(colMedians(rle)^2) # Variance of the median

-  

-  # Variable IQR RLE

-  RLE_IQR <- var(colIQRs(rle)) # Variance of the IQR

+  if (stratified_rle) {

+    biobatch = as.factor(paste(bio, batch, sep = "_"))

+    RLE_MED <- RLE_IQR <- 0

+    

+    for (cond in levels(biobatch)) {

+      is_cond = which(biobatch == cond)

+      cond_w = length(is_cond)

+      rle <- expr[, is_cond] - rowMedians(expr[, is_cond])

+      

+      # Non-Zero Median RLE

+      RLE_MED <-

+        RLE_MED + cond_w * mean(colMedians(rle) ^ 2) # Var of the med

+      

+      # Variable IQR RLE

+      RLE_IQR <-

+        RLE_IQR + cond_w * var(colIQRs(rle)) # Variance of the IQR

+      

+    }

+    

+    RLE_MED <- RLE_MED / length(biobatch)

+    RLE_IQR <- RLE_IQR / length(biobatch)

+    

+  } else{

+    rle <- expr - rowMedians(expr)

+    

+    # Non-Zero Median RLE

+    RLE_MED <- mean(colMedians(rle) ^ 2) # Variance of the median

+    

+    # Variable IQR RLE

+    RLE_IQR <- var(colIQRs(rle)) # Variance of the IQR

+  }

-  scores = c(BIO_SIL, BATCH_SIL, PAM_SIL,

-             EXP_QC_COR, EXP_UV_COR, EXP_WV_COR,

-             RLE_MED, RLE_IQR)

-  names(scores) = c("BIO_SIL", "BATCH_SIL", "PAM_SIL",

-                    "EXP_QC_COR", "EXP_UV_COR", "EXP_WV_COR",

-                    "RLE_MED", "RLE_IQR")

+  scores = c(BIO_SIL,

+             BATCH_SIL,

+             PAM_SIL,

+             EXP_QC_COR,

+             EXP_UV_COR,

+             EXP_WV_COR,

+             RLE_MED,

+             RLE_IQR)

+  names(scores) = c(

+    "BIO_SIL",

+    "BATCH_SIL",

+    "PAM_SIL",

+    "EXP_QC_COR",

+    "EXP_UV_COR",

+    "EXP_WV_COR",

+    "RLE_MED",

+    "RLE_IQR"

+  )

   return(scores)

 }

@@ -50,7 +50,13 @@

 #' @param stratified_pam logical. If TRUE then maximum ASW for PAM_SIL is

 #'   separately computed for each biological-cross-batch stratum (accepting

 #'   NAs), and a weighted average is returned as PAM_SIL.

-#'   

+#' @param stratified_cor logical. If TRUE then cor metrics are separately

+#'   computed for each biological-cross-batch stratum (accepts NAs), and

+#'   weighted averages are returned for EXP_QC_COR, EXP_UV_COR, & EXP_WV_COR.

+#'   Default FALSE.

+#' @param stratified_rle logical. If TRUE then rle metrics are separately 

+#'   computed for each biological-cross-batch stratum (accepts NAs), and

+#'   weighted averages are returned for RLE_MED & RLE_IQR. Default FALSE.

 #' @param verbose logical. If TRUE some messagges are printed.

 #' @param run logical. If FALSE the normalization and evaluation are not run,

 #'   but normalization parameters are returned in the output object for 

@@ -169,14 +175,16 @@ setMethod(

                         eval_proj = NULL,

                         eval_proj_args = NULL, eval_kclust=2:10,

                         verbose=FALSE, stratified_pam = FALSE,

+                        stratified_cor = FALSE,

+                        stratified_rle = FALSE,

                         return_norm = c("no", "in_memory", "hdf5"),

                         hdf5file,

                         bpparam=BiocParallel::bpparam()) {

-    

-    

-    rezero = (zero %in% c("preadjust","strong"))

-    fixzero = (zero %in% c("postadjust","strong"))

-    

+

+    zero <- match.arg(zero)

+    rezero <- (zero %in% c("preadjust","strong"))

+    fixzero <- (zero %in% c("postadjust","strong"))

+

     if(x@is_log) {

       stop("At the moment, scone is implemented only for non-log counts.")

     }

@@ -329,6 +337,18 @@ setMethod(

         stop("'eval_kclust' must be less than the number of samples.")

       }

+      if(stratified_rle) {

+        if(is.null(bio) & is.null(batch)){

+          stop("For stratified_rle, bio and/or batch must be specified")

+        }

+      }

+      

+      if(stratified_cor) {

+        if(is.null(bio) & is.null(batch)){

+          stop("For stratified_cor, bio and/or batch must be specified")

+        }

+      }

+      

       if(!is.null(eval_kclust) & stratified_pam) {

         if(is.null(bio) & is.null(batch)){

           stop("For stratified_pam, bio and/or batch must be specified")

@@ -539,7 +559,10 @@ setMethod(

                               bio = bio, batch = batch,

                               qc_factors = qc_pcs, 

                               uv_factors = uv_factors, wv_factors = wv_factors,

-                              is_log = TRUE, stratified_pam = stratified_pam)

+                              stratified_pam = stratified_pam,

+                              stratified_cor = stratified_cor,

+                              stratified_rle = stratified_rle,

+                              is_log = TRUE)

       } else {

         score <- NULL

       }

@@ -296,8 +296,8 @@ impute_null <- function(expression,impute_args) {

 #'   \item{p_nodrop}{ 1 - the probability P(drop|Y), useful as weights in 

 #'   weighted PCA} \item{expected_state}{ the expected

 #'   value E[Z] (1 = "on")} \item{loglik}{ the log-likelihood} 

-#'   \item{convergence}{for all genes, 0 if the algorithm converged and 1 if maxiter was 

-#'   reached} }

+#'   \item{convergence}{for all genes, 0 if the algorithm converged and

+#'   1 if maxiter was reached} }

 #'   

 #' @examples

 #' mat <- matrix(rpois(1000, lambda = 3), ncol=10)

@@ -1,3 +1,10 @@

+Changes in version 1.1.3 (2018-10-19)

+========================

+

+* Modified COR scores to reflect R^2 for regression on UV/WV/QC rather than max correlation.

+* Bug fix for shiny report.

+* Updated vignette.

+

 Changes in version 1.1.2 (2017-7-3)

 ========================

@@ -2,7 +2,8 @@

 % Please edit documentation in R/SCONE_DEFAULTS.R

 \name{DESEQ_FN}

 \alias{DESEQ_FN}

-\title{Relative log-expression (RLE; DESeq) scaling normalization wrapper function}

+\title{Relative log-expression (RLE; DESeq) scaling normalization wrapper

+function}

 \usage{

 DESEQ_FN(ei)

 }

@@ -13,7 +14,8 @@ DESEQ_FN(ei)

 RLE normalized matrix.

 }

 \description{

-Relative log-expression (RLE; DESeq) scaling normalization wrapper function

+Relative log-expression (RLE; DESeq) scaling normalization wrapper

+function

 }

 \details{

 SCONE scaling wrapper for \code{\link[edgeR]{calcNormFactors}}).

@@ -40,8 +40,8 @@ a list with the following elements: \itemize{ \item{W}{ coefficients

   \item{p_nodrop}{ 1 - the probability P(drop|Y), useful as weights in 

   weighted PCA} \item{expected_state}{ the expected

   value E[Z] (1 = "on")} \item{loglik}{ the log-likelihood} 

-  \item{convergence}{for all genes, 0 if the algorithm converged and 1 if maxiter was 

-  reached} }

+  \item{convergence}{for all genes, 0 if the algorithm converged and

+  1 if maxiter was reached} }

 }

 \description{

 This function implements Newton's method for solving zero of 

@@ -15,8 +15,8 @@ scone(x, ...)

   adjust_batch = c("no", "yes", "force"), run = TRUE, evaluate = TRUE,

   eval_pcs = 3, eval_proj = NULL, eval_proj_args = NULL,

   eval_kclust = 2:10, verbose = FALSE, stratified_pam = FALSE,

-  return_norm = c("no", "in_memory", "hdf5"), hdf5file,

-  bpparam = BiocParallel::bpparam())

+  stratified_cor = FALSE, stratified_rle = FALSE, return_norm = c("no",

+  "in_memory", "hdf5"), hdf5file, bpparam = BiocParallel::bpparam())

 }

 \arguments{

 \item{x}{a \code{\link{SconeExperiment}} object.}

@@ -77,6 +77,15 @@ If NULL, tightness will be returned NA.}

 separately computed for each biological-cross-batch stratum (accepting

 NAs), and a weighted average is returned as PAM_SIL.}

+\item{stratified_cor}{logical. If TRUE then cor metrics are separately

+computed for each biological-cross-batch stratum (accepts NAs), and

+weighted averages are returned for EXP_QC_COR, EXP_UV_COR, & EXP_WV_COR.

+Default FALSE.}

+

+\item{stratified_rle}{logical. If TRUE then rle metrics are separately 

+computed for each biological-cross-batch stratum (accepts NAs), and

+weighted averages are returned for RLE_MED & RLE_IQR. Default FALSE.}

+

 \item{return_norm}{character. If "no" the normalized values will not be

 returned with the output object. This will create a much smaller object

 and may be useful for large datasets and/or when many combinations are

@@ -7,10 +7,10 @@

 score_matrix(expr, eval_pcs = 3, eval_proj = NULL, eval_proj_args = NULL,

   eval_kclust = NULL, bio = NULL, batch = NULL, qc_factors = NULL,

   uv_factors = NULL, wv_factors = NULL, is_log = FALSE,

-  stratified_pam = FALSE)

+  stratified_pam = FALSE, stratified_cor = FALSE, stratified_rle = FALSE)

 }

 \arguments{

-\item{expr}{matrix. The expression data matrix (genes in rows, cells in 

+\item{expr}{matrix. The expression data matrix (genes in rows, cells in

 columns).}

 \item{eval_pcs}{numeric. The number of principal components to use for

@@ -19,11 +19,11 @@ evaluation (Default 3). Ignored if !is.null(eval_proj).}

 \item{eval_proj}{function. Projection function for evaluation (see Details).

 If NULL, PCA is used for projection}

-\item{eval_proj_args}{list. List of arguments passed to projection function 

+\item{eval_proj_args}{list. List of arguments passed to projection function

 as eval_proj_args (see Details).}

-\item{eval_kclust}{numeric. The number of clusters (> 1) to be used for pam 

-tightness (PAM_SIL) evaluation. If an array of integers, largest average 

+\item{eval_kclust}{numeric. The number of clusters (> 1) to be used for pam

+tightness (PAM_SIL) evaluation. If an array of integers, largest average

 silhouette width (tightness) will be reported in PAM_SIL. If NULL, PAM_SIL

 will be returned NA.}

@@ -33,10 +33,10 @@ NA is allowed. If NULL, condition ASW, BIO_SIL, will be returned NA.}

 \item{batch}{factor. A known batch variable (variation to be removed), NA is

 allowed. If NULL, batch ASW, BATCH_SIL, will be returned NA.}

-\item{qc_factors}{Factors of unwanted variation derived from quality 

+\item{qc_factors}{Factors of unwanted variation derived from quality

 metrics. If NULL, qc correlations, EXP_QC_COR, will be returned NA.}

-\item{uv_factors}{Factors of unwanted variation derived from negative 

+\item{uv_factors}{Factors of unwanted variation derived from negative

 control genes (evaluation set). If NULL, uv correlations, EXP_UV_COR,

 will be returned NA.}

@@ -44,36 +44,47 @@ will be returned NA.}

 control genes (evaluation set). If NULL, wv correlations, EXP_WV_COR,

 will be returned NA.}

-\item{is_log}{logical. If TRUE the expr matrix is already logged and log 

+\item{is_log}{logical. If TRUE the expr matrix is already logged and log

 transformation will not be carried out prior to projection. Default FALSE.}

-\item{stratified_pam}{logical. If TRUE then maximum ASW is separately 

-computed for each biological-cross-batch stratum (accepts NAs), and a 

+\item{stratified_pam}{logical. If TRUE then maximum ASW is separately

+computed for each biological-cross-batch stratum (accepts NAs), and a

 weighted average silhouette width is returned as PAM_SIL. Default FALSE.}

+

+\item{stratified_cor}{logical. If TRUE then cor metrics are separately

+computed for each biological-cross-batch stratum (accepts NAs), and

+weighted averages are returned for EXP_QC_COR, EXP_UV_COR, & EXP_WV_COR.

+Default FALSE.}

+

+\item{stratified_rle}{logical. If TRUE then rle metrics are separately

+computed for each biological-cross-batch stratum (accepts NAs), and

+weighted averages are returned for RLE_MED & RLE_IQR. Default FALSE.}

 }

 \value{

 A list with the following metrics: \itemize{ \item{BIO_SIL}{ Average

-  silhouette width by biological condition.} \item{BATCH_SIL}{ Average 

-  silhouette width by batch condition.} \item{PAM_SIL}{ Maximum average 

-  silhouette width from PAM clustering (see stratified_pam argument).} 

-  \item{EXP_QC_COR}{ Maximum squared spearman correlation between expression

-  pcs and quality factors.} \item{EXP_UV_COR}{ Maximum squared spearman