@@ -24,11 +24,13 @@ Imports:

   fpc,

   gplots,

   grDevices,

+  hexbin,

   limma,

   MASS,

   matrixStats,

   mixtools,

   grDevices,

+  RColorBrewer,

   boot,

   shiny,

   miniUI,

@@ -21,12 +21,14 @@ export(lm_adjust)

 export(make_design)

 export(metric_sample_filter)

 export(scone)

+export(scone_wrapper)

 export(score_matrix)

 import(BiocParallel)

 import(gplots)

 importFrom(DESeq,estimateSizeFactorsForMatrix)

 importFrom(EDASeq,betweenLaneNormalization)

 importFrom(MASS,glm.nb)

+importFrom(RColorBrewer,brewer.pal)

 importFrom(RUVSeq,RUVg)

 importFrom(aroma.light,normalizeQuantileRank.matrix)

 importFrom(boot,inv.logit)

@@ -44,6 +46,8 @@ importFrom(graphics,hist)

 importFrom(graphics,par)

 importFrom(graphics,plot)

 importFrom(graphics,text)

+importFrom(hexbin,hexbin)

+importFrom(hexbin,plot)

 importFrom(limma,lmFit)

 importFrom(matrixStats,colIQRs)

 importFrom(matrixStats,colMedians)

deleted file mode 100644

@@ -1,420 +0,0 @@

-library(Biobase) #for AnnotatedDataFrame if scone is not loaded

-library(scone)

-library(RColorBrewer)

-library(boot) #for logit

-

-#A wrapper for filtering and running scone

-#currently supports only one-variable bio and batch

-scone_something = function(selected_eSet, housekeeping_list,

-                           batch, bio = NULL, de = NULL,

-                           adjust_bio = "no", adjust_batch = "yes", rezero = TRUE,

-                           out_dir = getwd(), stratified_pam = TRUE,

-                           selected_expression_in_eSet = "counts_table",

-                           results_image_filename = "sconeResults.Rda",

-                           plots_filename = "scone_plots.pdf") {

-

-printf <- function(...) cat(sprintf(...))

-set.seed(112233) ## reproducibility

-if (!file.exists(out_dir)) {

-  dir.create(out_dir, recursive=TRUE)

-}

-pdf(file.path(out_dir, plots_filename), onefile=TRUE)

-

-#keep only protein coding genes

-keepGenes = featureData(selected_eSet)$Transcript_Type == "protein_coding"

-selected_eSet = selected_eSet[keepGenes,]

-printf("Removed %d non-protein coding genes, %d protein-coding genes left\n", sum(!keepGenes), sum(keepGenes))

-

-

-#don't use exprs because we have several units (e.g., counts, tpm, fpkm) for each eSet

-#instead, dynamically access the assayData with the user's argument

-selected_unit_exprs = get(selected_expression_in_eSet, envir = assayData(selected_eSet))

-selected_unit_exprs[is.na(selected_unit_exprs)] = 0

-E = aggregate(selected_unit_exprs, by=list(featureData(selected_eSet)$Gene_Symbol),

-              FUN=max, na.rm=TRUE)

-rownames(E) = E[,1]

-E = E[, -1]

-E = data.matrix(E)

-E[E == -Inf] = 0 #cases in which the max had no non-NA option - should not happen because we turned NAs into 0's already

-stopifnot(all(is.finite(E)))

-printf("collapsed %d transcripts into %d genes by gene symbols\n", nrow(selected_unit_exprs), nrow(E))

-

-#now to remove duplicated features from the featuresData (we know that the same gene

-#symbol will have the same type, so no harm in just removing duplicate rows) and then

-#reordering it according to the new expression matrix

-featureDat = featureData(selected_eSet)

-featureDat = featureDat[!duplicated(featureDat$Gene_Symbol)]

-featureDat = featureDat[rownames(E), ]

-

-# Expression set object to make filtering easier

-chosen_units_eSet = new("ExpressionSet", exprs = E, featureData = featureDat, 

-                        protocolData = protocolData(selected_eSet), phenoData = phenoData(selected_eSet))

-

-# Remove failed samples:

-is.failed = apply(is.na(exprs(chosen_units_eSet)) | (exprs(chosen_units_eSet) == 0), 2, all) |

-  apply(is.na(pData(protocolData(chosen_units_eSet))), 1, all)

-chosen_units_eSet = chosen_units_eSet[,!is.failed ]

-printf("Removed %d failed samples\n", sum(is.failed))

-

-

-## ----- Pre-Filtering of Transcripts: ----

-# Select only coding transcripts and TCR segments

-is.expressed.sc = rowMeans(exprs(chosen_units_eSet)) > 0

-chosen_units_eSet = chosen_units_eSet[is.expressed.sc, ]

-printf("Removed %d undetected genes, %d left\n",sum(!is.expressed.sc), sum(is.expressed.sc))

-

-

-# Params for gene filtering

-# turns out that the filtering of the Th17 dataset was pretty robust to the choice of the quantile

-thresh_fail = quantile(exprs(chosen_units_eSet)[exprs(chosen_units_eSet) > 0], 0.2) #10

-num_fail = 10

-

-# Initial Gene Filtering

-init.gf.vec = rowSums(exprs(chosen_units_eSet) > thresh_fail) > num_fail

-chosen_units_eSet = chosen_units_eSet[init.gf.vec, ]

-printf("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))

-

-

-

-

-

-

-#The counts data frame contains feature-level read counts from tophat alignments of 96 single-cell libraries to the mm10 reference genome [@trapnell2009]. qc contains library alignment metrics obtained from Picard.

-#de and hk are positive and negative control gene sets derived from population studies. batch and bio are factors labeling batch and time point respectively. Consider the joint distribution of these factors:

-counts = exprs(chosen_units_eSet)

-qc = pData(protocolData(chosen_units_eSet))

-

-#prevent casing problems

-rownames(counts) = sapply(rownames(counts), toupper)

-housekeeping_genes = sapply(as.vector(as.matrix(read.table(housekeeping_list, stringsAsFactors = FALSE))), toupper)

-hk = intersect(housekeeping_genes, rownames(counts))

-

-

-

-#####################################

-#start to do real work

-#####################################

-unnecessaryQC_metrics = apply(is.na(qc), 2, any) | apply(qc, 2, sd) < 1e-5

-printf("remove %d QC metrics (%s) that either contained NAs or were constant across cells\n", sum(unnecessaryQC_metrics), paste(colnames(qc)[unnecessaryQC_metrics], collapse = "; "))

-qc = qc[, !unnecessaryQC_metrics]

-

-colnames(qc)

-

-## Joint distribution of batches and biological conditions (time after induction)

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

-  table(batch,bio)

-}

-

-# Color scheme

-cc <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"),brewer.pal(9,"Set3"))

-

-# Barplot of read proportion mapping to the genome

-barplot(qc$RALIGN[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Percentage of mapped reads, colored by batch")

-legend("bottomleft", legend=levels(batch), fill=cc,cex=0.4)

-

-# Barplot of total read number

-barplot(qc$NREADS[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Total number of reads, colored by batch")

-legend("topright", legend=levels(batch), fill=cc, cex=0.4)

-

-qpc = prcomp(qc, center = TRUE, scale = TRUE)

-barplot(cumsum(qpc$sdev^2)/sum(qpc$sdev^2), border="gray", xlab="PC", ylab="proportion of variance", main="Quality PCA")

-

-barplot(qpc$x[,1][order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="Quality PC1, colored by batch")

-legend("bottomleft", legend=levels(batch), fill=cc, cex=0.8)

-

-

-# Mean log10(x+1) expression

-mu_obs = rowMeans(log10(counts[hk,]+1)) #0 if hk == NULL

-

-# Drop-outs

-drop_outs = counts[hk,] == 0

-

-# Logistic Regression Model of Failure

-logistic_coef = matrix(0,ncol(drop_outs),2)

-for (si in seq_len(ncol(drop_outs))){

-  fit = suppressWarnings(glm(cbind(drop_outs[,si],1 - drop_outs[,si]) ~ mu_obs,family=binomial(logit)))

-  if(fit$converged){

-    logistic_coef[si,] = fit$coefficients

-  } else {

-    logistic_coef[si,1] = logit(drop_rate[si])

-  }

-}

-

-

-par(mfrow=c(1,2))

-

-# Plot Failure Curves and Calculate AUC

-plot(NULL, main = "False Negative Rate Curves",ylim = c(0,1),xlim = c(0,6), ylab = "Failure Probability", xlab = "Mean log10 Expression")

-x = (0:60)/10

-AUC = NULL

-for(si in 1:ncol(counts)){

-  y = 1/(exp(-logistic_coef[si,1] - logistic_coef[si,2] * x) + 1)

-  AUC[si] = sum(y)/10

-  lines(x, 1/(exp(-logistic_coef[si,1] - logistic_coef[si,2] * x) + 1), type = 'l', lwd = 2, col = cc[batch][si])

-}

-

-# Barplot of FNR AUC

-barplot(AUC[order(batch)], col=cc[batch][order(batch)], border=cc[batch][order(batch)], main="FNR AUC, colored by batch")

-legend("topleft", legend=levels(batch), fill=cc, cex=0.4)

-

-mfilt_report <- metric_sample_filter(expr = counts,

-                                     nreads = qc$NREADS,ralign = qc$RALIGN,

-                                     suff_nreads = 10^5,

-                                     suff_ralign = 90,

-                                     pos_controls = rownames(counts) %in% hk,

-                                     zcut = 3,mixture = FALSE, plot = TRUE)

-

-

-hist(qc$RALIGN, breaks = 0:100)

-# Hard threshold

-abline(v = 15, col = "yellow", lwd = 2) 

-# 3 (zcut) standard deviations below the mean ralign value

-abline(v = mean(qc$RALIGN) - 3*sd(qc$RALIGN), col = "green", lwd = 2) 

-# 3 (zcut) MADs below the median ralign value

-abline(v = median(qc$RALIGN) - 3*mad(qc$RALIGN), col = "red", lwd = 2)

-# Sufficient threshold

-abline(v = 90, col = "grey", lwd = 2)

-# Final threshold is the minimum of 1) the sufficient threshold and 2) the max of all others

-thresh = min(90,max(c(15,mean(qc$RALIGN) - 3*sd(qc$RALIGN),median(qc$RALIGN) - 3*mad(qc$RALIGN))))

-abline(v = thresh, col = "blue", lwd = 2, lty = 2)

-

-legend("topleft",legend = c("Hard","SD","MAD","Sufficient","Final"),lwd = 2, col = c("yellow","green","red","grey","blue"),lty = c(1,1,1,1,2), cex = .5)

-

-

-# Which thresholds are missing? (breadth)

-is_na_filt = unlist(lapply(is.na(mfilt_report),any)) 

-

-# Identify samples failing any threshold

-m_sampfilter = !apply(simplify2array(mfilt_report[!is_na_filt]),1,any)

-

-# Filter Samples

-fcounts = counts[,m_sampfilter]

-fqc = qc[m_sampfilter,]

-fbatch = batch[m_sampfilter]

-fbio = bio[m_sampfilter] #if bio==NULL then we get from this line fbio=NULL

-

-# # Simple gene filter --> redundant with what I have above

-# filterCount <- function(counts, nRead=5, nCell=5){

-#   filter <- apply(counts, 1, function(x) length(x[x>=nRead])>=nCell)

-#   return(filter)

-# }

-# genefilter <- filterCount(fcounts)

-# fcounts = fcounts[genefilter,]

-

-# skip filtering genes again

-fcounts = counts

-fhk = hk[hk %in% rownames(fcounts)]

-fde = de[de %in% rownames(fcounts)] #if de==NULL then we get from this line fde=NULL

-

-

-#FNR

-#Michael: I see you recommend bulk_model=true even though the default is false? Why?

-#I'm not clear on what this parameter does - what other gene features are possible if gFeatM is not specified?

-fnr_out = estimate_ziber(x = fcounts, bulk_model = TRUE,

-                         pos_controls = rownames(fcounts) %in% fhk,

-                         verbose = TRUE, maxiter = 1000)

-

-

-SUM_FN = function (ei) 

-{

-  sums = colSums(ei)

-  eo = t(t(ei)*mean(sums)/sums)

-  return(eo)

-}

-

-MED_FN_POS = function (ei) 

-{

-  zei = ei

-  zei[ei == 0] = NA

-  q = apply(zei, 2, quantile, 0.5, na.rm = TRUE)

-  zeo = t(t(zei)/q) * mean(q)

-  eo = zeo

-  eo[ei == 0] = 0

-  return(eo)

-}

-

-EFF_FN = function (ei) 

-{

-  sums = colSums(ei > 0)

-  eo = t(t(ei)*sums/mean(sums))

-  return(eo)

-}

-

-imputation = list(none=impute_null,expect=impute_expectation)

-impute_args = list(p_nodrop = fnr_out$p_nodrop, mu = exp(fnr_out$Alpha[1,]))

-

-scaling = list(none=identity,

-             sum = SUM_FN,

-             eff = EFF_FN,

-             tmm = TMM_FN,

-             uq = UQ_FN,

-             uqp = UQ_FN_POS,

-             fq = FQT_FN,

-             fqp = FQ_FN_POS,

-             deseq=DESEQ_FN,

-             deseqp=DESEQ_FN_POS)

-

-

-

-# Generate Params (RUN = FALSE)

-#Michael: in your example qc=ppq which is ppq = scale(qc,center = TRUE,scale = TRUE)

-#but I see in scone's code that it expects QC matrix and does the prcomp itself

-params <- scone(expr = as.matrix(fcounts), scaling = scaling,

-                imputation = imputation, impute_args = impute_args,

-                ruv_negcon = fhk, k_ruv = 3,

-                qc = as.matrix(fqc), k_qc = 3,

-                bio = fbio, adjust_bio = adjust_bio,

-                batch = fbatch, adjust_batch = adjust_batch,

-                run = FALSE)

-head(params)

-

-apply(params,2,unique)

-

-#exclide batch adjustment and no adjustment by bio and similar not-making-sense scenarios

-is_screened = ((params$imputation_method == "none") & (params$scaling_method %in% c("deseq","uq","tmm"))) | 

-  ((params$imputation_method == "expect") & (params$scaling_method %in% c("deseqp","uqp","fqp","eff"))) |

-  ((params$adjust_biology == "bio") & (params$adjust_batch != "batch"))

-

-params = params[!is_screened,]

-

-##NO MORE PLOTS!

-dev.off()

-

-# Generate Scores and Ranking

-tic = proc.time()

-res = scone(expr = as.matrix(fcounts), scaling = scaling,

-            imputation = imputation, impute_args = impute_args,

-            ruv_negcon = fhk, k_ruv = 3,

-            qc = as.matrix(fqc), k_qc = 3,

-            bio = fbio, adjust_bio = adjust_bio,

-            batch = fbatch, adjust_batch = adjust_batch,

-            run = TRUE, params = params, verbose = TRUE,

-            eval_poscon = fde, eval_kclust = 2:3, stratified_pam = stratified_pam,

-            rezero = rezero, return_norm = "in_memory")

-toc = proc.time()

-

-printf("time elapsed:\n")

-print(toc-tic)

-

-save.image(file.path(out_dir, results_image_filename))

-

-

-printf("Scone finished successfully and image saved")

-

-}

-

-

-#code to run this method and test it

-if(FALSE) {

-  rm(list=ls())

-  #library(bioc2016singlecell)

-  set.seed(112233) ## for reproducibility

-  printf <- function(...) cat(sprintf(...))

-  

-  library(Biobase) #for AnnotatedDataFrame if scone is not loaded

-  library(scone)

-  library(RColorBrewer)

-  

-  BiocParallel::register(BiocParallel::SerialParam())

-  

-  setwd("/data/yosef2/Published_Data/TH17/code")

-  source("scone_something.R")

-  

-  load("/data/yosef2/Published_Data/TH17/code/combined_th17_metadata.Rdata")

-  housekeeping_list = "/data/yosef/CD8_effector_diff/src/SummaryPipeline/house_keeping_mouse_TitleCase.txt"

-  out_dir = "/data/yosef2/Published_Data/TH17processed_aw20160712/scone/only_T16"

-  

-  #only T16 from the non-GFP mouse

-  keep = combinedMetadata$michael_condition == "TGFB1_IL6-48h" &

-    !combinedMetadata$is_population

-  

-  combinedMetadata = combinedMetadata[keep,]

-  

-  #tophaht + featureCounts

-  collect_dir="/data/yosef2/Published_Data/TH17/processed_aw20160712/collect"

-  load(file.path(collect_dir, "collectedRNASeqStudy.RData"))

-  selected_eSet = collectedRNASeqStudy$cuff_eSet

-  #featureNames(selected_eSet) = toupper(featureNames(selected_eSet))

-  

-  selected_eSet = selected_eSet[, combinedMetadata$SRR]

-  batch = combinedMetadata$michael_batch

-  scone_something(selected_eSet, housekeeping_list = housekeeping_list,

-                  selected_expression_in_eSet = "counts_table",

-                  batch = batch, bio = NULL, de = NULL,

-                  adjust_bio = "no", adjust_batch = "yes", rezero = TRUE,

-                  out_dir = getwd(), stratified_pam = TRUE,

-                  results_image_filename = "sconeResults_tophatFeatureCounts.Rda",

-                  plots_filename = "scone_plot_tophatFeatureCounts.pdf")

-  

-  stop("finished")

-  

-}

-

-

-#code to inspect scone results

-if(FALSE) {

-  rm(list=ls())

-  set.seed(4455)

-  out_dir = "/home/eecs/allonwag/data2/TFH/processed_20160628/scone"

-  load(file.path(out_dir, "sconeResults.Rda"))

-  

-  

-  names(res)

-  

-  if(is.null(bio)) {

-    #no positive control genes (de) so no evaluation of positive controls

-    BIO_DEPENDENT_EVALUATION_METRICS = c("BIO_SIL", "EXP_WV_COR")

-    res$scores = res$scores[, !(colnames(res$scores) %in% BIO_DEPENDENT_EVALUATION_METRICS)]

-    res$metrics = res$metrics[, !(colnames(res$metrics) %in% BIO_DEPENDENT_EVALUATION_METRICS)]

-  }

-  

-  head(res$scores)

-  

-  pc_obj = prcomp(res$scores[,-ncol(res$scores)],center = TRUE,scale = FALSE)

-  bp_obj = biplot_colored(pc_obj,y = -res$scores[,ncol(res$scores)],expand = .6)

-  

-  #library(scone)

-  #if I don't have this library installed

-  source("biplot.R")

-  bp_obj = biplot_colored(pc_obj,y = -res$scores[,ncol(res$scores)],expand = .6)

-  

-  points(bp_obj[grepl(",bio,batch",rownames(bp_obj)),], pch = 1, col = "red", cex = 1)

-  points(bp_obj[grepl(",bio,batch",rownames(bp_obj)),], pch = 1, col = "red", cex = 1.5)

-  

-  

-  bp_obj = biplot_colored(pc_obj,y = -res$scores[,ncol(res$scores)],expand = .6)

-  

-  points(bp_obj[grepl(",no_bio,batch",rownames(bp_obj)),], pch = 1, col = "red", cex = 1)

-  points(bp_obj[grepl(",no_bio,batch",rownames(bp_obj)),], pch = 1, col = "red", cex = 1.5)

-  

-  #my addition to scone's tutorial

-  bp_obj = biplot_colored(pc_obj,y = -res$scores[,ncol(res$scores)],expand = .6)

-  

-  points(bp_obj[grepl(",fq,",rownames(bp_obj)),], pch = 1, col = "green", cex = 1.5)

-  points(bp_obj[grepl(",tmm,",rownames(bp_obj)),], pch = 1, col = "cyan", cex = 1.5)

-  points(bp_obj[grepl(",deseq,",rownames(bp_obj)),], pch = 1, col = "magenta", cex = 1.5)

-  points(bp_obj[grepl(",uqp,",rownames(bp_obj)),], pch = 1, col = "black", cex = 1.5)

-  #end my addition

-  

-  bp_obj = biplot_colored(pc_obj,y = -res$scores[,ncol(res$scores)],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) == rownames(params)[1],]), pch = 1, col = "blue", cex = 1)

-  points(t(bp_obj[rownames(bp_obj) == rownames(params)[1],]), pch = 1, col = "blue", cex = 1.5)

-  

-  arrows(bp_obj[rownames(bp_obj) == rownames(params)[1],][1],

-         bp_obj[rownames(bp_obj) == rownames(params)[1],][2],

-         bp_obj[1,][1],

-         bp_obj[1,][2],

-         lty = 2, lwd = 2)

-  

-  printf("(Trumpets): the selected method is %s", rownames(res$scores)[1])

-  normalized_expression = res$normalized_data[[1]]

-  selected_normalization_method = names(res$normalized_data)[1]

-  save(normalized_expression, selected_normalization_method, file=file.path(out_dir, "sconeNormalizedExpression.Rda"))

-  

-  sessionInfo()

-  

-}

new file mode 100644

@@ -0,0 +1,236 @@

+#' Wrapper for Running Essential SCONE Modules

+#'

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

+#' @param qc matrix. 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.

+#'   Default "0".

+#' @param out_dir character. Output directory. 

+#'   Default getwd().

+#' @param seed numeric. Random seed. 

+#'   Default 112233.

+#'   

+#' @param filt_genes logical. Should genes be filtered post-sample filtering?

+#'   Default TRUE.

+#'   

+#' @param fnr_maxiter numeric. Maximum number of iterations in EM estimation of expression posteriors.

+#'    Dafault 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("sum", "deseq", "tmm", "uq", "fq"). See details.

+#' @param norm_rezero logical. Restore imputed zeroes to zero following the Scaling Step?

+#'   Default TRUE.

+#' @param norm_k_max numeric. Maximum 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.2. 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_kclust numeric. The number of clusters (> 1) to be used for pam

+#'   tightness evaluation. If an array of integers, largest average silhouette

+#'   width (tightness) will be reported. 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.

+#'

+#' @export

+#' @importFrom RColorBrewer brewer.pal

+#' @importFrom boot inv.logit

+#' @importFrom hexbin plot hexbin

+#'

+#' @details "ADD DESCRIPTION"

+#'

+#' @return Directory structure "ADD DESCRIPTION"

+#'

+

+scone_wrapper <- function(expr, qc,

+                  bio=NULL, batch=NULL, negcon=NULL,

+                  verbose=c("0","1","2"), out_dir=getwd(), seed = 112233,

+                  

+                  filt_genes=TRUE,

+                  

+                  fnr_maxiter=1000,

+                  

+                  norm_impute=c("yes", "no", "force"),

+                  norm_scaling=c("sum", "deseq", "tmm", "uq", "fq"),

+                  norm_rezero=TRUE, 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.2,

+                  eval_poscon=NULL, eval_kclust = 2:10, eval_stratified_pam = TRUE) {

+

+  verbose = match.arg(verbose)

+  verbose = as.numeric(verbose)

+  

+  printf <- function(...) cat(sprintf(...))

+  set.seed(seed)

+  

+  if(!is.matrix(expr)){stop("expr must be a matrix")}

+  

+  # 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)

+  }

+  

+  ## ------ Data Filtering Module ------

+  

+  # 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("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("Data Filtering Module: %d negative control genes to be used for sample filtering\n", length(small_negcon))}

+  

+  # Metric-based Filtering and Report

+  pdf(paste0(misc_dir,"/filter_report.pdf"))

+  mfilt = metric_sample_filter(expr, mixture = TRUE, plot = TRUE, hist_breaks = 100,

+                               zcut = 4, ###Make argument###

+                               pos_controls = rownames(ct) %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("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,]

+  batch = droplevels(batch[mfilt])

+  bio = droplevels(bio[mfilt])

+

+  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(verbose > 0){printf("Data Filtering Module: Kept only %d genes expressed in more than %.2f units in more than %d cells, excluded %d genes\n", sum(final.gf.vec), thresh_fail, num_fail, sum(!final.gf.vec))}

+  

+    # Implement gene filter

+    expr = expr[final.gf.vec,]

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

+      negcon = negcon[negcon %in% rownames(expr)]

+      if(verbose > 0){printf("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("Data Filtering Module: Kept only %d positive control genes\n", length(eval_poscon))}

+    }

+  }

+  

+  if(verbose > 0){printf("Data Filtering Module: COMPLETE\n\n")}

+  ## ------ False-Negative Rate Inference Module ------

+  

+  if(verbose > 0){printf("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% hk,

+                           verbose = (verbose > 1), 

+                           maxiter = fnr_maxiter)

+  if(fnr_out$convergence == 1){printf("False-Negative Rate Inference Module: WARNING - EM did not converge. Try increasing fnr_max_iter\n")}

+  if(verbose > 0){printf("False-Negative Rate Inference Module: Producing report...\n")}

+  

+  # FNR Report

+  cc <- c(brewer.pal(9, "Set1"), brewer.pal(8, "Set2"),brewer.pal(9,"Set3"))

+  

+  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)))

+    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("False-Negative Rate Inference Module: COMPLETE\n\n")}

+  

+  dev.off()

+  

+  ui <- fluidPage(

+    h1("Example app"),

+    sidebarLayout(

+      sidebarPanel(

+        actionButton("qtBtn", "Quit")

+      ),

+      mainPanel(

+      )

+    )

+  )

+

+  server <- function(input, output, session) {

+    

+    observeEvent(input$qtBtn,{

+      if(input$qtBtn > 0){

+        stopApp("hello world")

+      }

+    })

+  }

+  

+  print(runApp(list(ui = ui, server = server)))

+

+  sessionInfo()

+  

+}

\ No newline at end of file

new file mode 100644

@@ -0,0 +1,107 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/scone_wrap.R

+\name{scone_wrapper}

+\alias{scone_wrapper}

+\title{Wrapper for Running Essential SCONE Modules}

+\usage{

+scone_wrapper(expr, qc, bio = NULL, batch = NULL, negcon = NULL,

+  verbose = c("0", "1", "2"), out_dir = getwd(), seed = 112233,

+  filt_genes = TRUE, fnr_maxiter = 1000, norm_impute = c("yes", "no",

+  "force"), norm_scaling = c("sum", "deseq", "tmm", "uq", "fq"),

+  norm_rezero = TRUE, 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.2,

+  eval_poscon = NULL, eval_kclust = 2:10, eval_stratified_pam = TRUE)

+}

+\arguments{

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

+

+\item{qc}{matrix. The quality control (QC) matrix (cells in rows, metrics in columns)

+to be used for filtering, normalization, and evaluation.}

+

+\item{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.}

+

+\item{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.}

+

+\item{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.}

+

+\item{verbose}{character. Verbosity level.

+Default "0".}

+

+\item{out_dir}{character. Output directory. 

+Default getwd().}

+

+\item{seed}{numeric. Random seed. 

+Default 112233.}

+

+\item{filt_genes}{logical. Should genes be filtered post-sample filtering?

+Default TRUE.}

+

+\item{fnr_maxiter}{numeric. Maximum number of iterations in EM estimation of expression posteriors.

+Dafault 1000.}

+

+\item{norm_impute}{character. Should imputation be included in the comparison? 

+If 'force', only imputed normalizations will be run.

+Default "yes."}

+

+\item{norm_scaling}{character. Scaling options to be included in the Scaling Step. 

+Default c("sum", "deseq", "tmm", "uq", "fq"). See details.}

+

+\item{norm_rezero}{logical. Restore imputed zeroes to zero following the Scaling Step?

+Default TRUE.}

+

+\item{norm_k_max}{numeric. Maximum number (norm_k_max) of factors of unwanted variation modeled

+in the Adjustment Step.

+Default NULL.}

+

+\item{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.}

+

+\item{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."}

+

+\item{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."}

+

+\item{eval_dim}{numeric. The number of principal components to use for evaluation.

+Default NULL.}

+

+\item{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.2. Ignored if eval_dim is not NULL.}

+

+\item{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.}

+

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

+tightness evaluation. If an array of integers, largest average silhouette

+width (tightness) will be reported. If NULL, tightness will be returned NA.}

+

+\item{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.}

+}

+\value{

+Directory structure "ADD DESCRIPTION"

+}

+\description{

+Wrapper for Running Essential SCONE Modules

+}

+\details{

+"ADD DESCRIPTION"

+}

+