diff --git a/R/calZscore.R b/R/calZscore.R
index e17c093..4074503 100644
--- a/R/calZscore.R
+++ b/R/calZscore.R
@@ -20,12 +20,18 @@ calZscore <- function(...) UseMethod("calZscore")
 
 calZscore.default <- function(x, mad = T) {
   if (mad) {
-    dev = mad(x, na.rm = T)
+    dev <- mad(x, na.rm = T)
   } else {
-    dev = sd(x, na.rm = T)
+    dev <- sd(x, na.rm = T)
   }
-  x = round((x-median(x, na.rm = T))/dev, 3)
-  x
+  # Avoid division by zero
+  if (dev  == 0) {
+    warning("Standard deviation is zero. Returning original values.")
+    return(x)
+  }
+
+  x <- round((x-median(x, na.rm = T))/dev, 3)
+  return(x)
 }
 
 #' Calculate z-scores of QC metrics of a SampleDataset
@@ -52,6 +58,7 @@ calZscore.sampleDataset <- function (
   object$zscore = .zscoreColNames(qcMetrics)
   object$df[object$zscore] = NA
 
+
   if (is.null(strat)) {
     # calculate z-score across all qcMetrics
     object$df[, object$zscore] = apply(object$df[, qcMetrics], 2, calZscore)
@@ -75,15 +82,26 @@ calZscore.sampleDataset <- function (
   return(object)
 }
 
-#' Calculate maximum vcf z-score and update the object with a maxVcfzscore
-#' column
+#' Calculate maximum VCF z-score and update the object with a maxVcfZscore column
+#'
+#' @param object a sampleDataset object
+#' @param qcMetrics a vector of names of QC metrics
+#' @return updated sampleDataset object with a new maxVcfZscore column
+#'
+#' @export
 .calMaxVcfZ <- function(object, qcMetrics) {
   vcfQcMetrInUse = .zscoreColNames(intersect(object$vcfQcMetr, qcMetrics))
   object$df$maxVcfzscore = apply(abs(object$df[vcfQcMetrInUse]), 1, max)
   return(object)
 }
 
+
 #' Get z-score column names
+#'
+#' @param qcMetrics a vector of QC metric names
+#' @return a vector of z-score column names
+#'
+#' @export
 .zscoreColNames <- function (qcMetrics) {
   return(paste(qcMetrics, 'Zscore', sep = ''))
 }
diff --git a/R/examples.R b/R/examples.R
deleted file mode 100644
index c969b0c..0000000
--- a/R/examples.R
+++ /dev/null
@@ -1,12 +0,0 @@
-#' Load example data
-load_examples <- function () {
-  bamQcMetr = read.csv('bamQcMetr.tsv', sep = '\t')
-  vcfQcMetr = read.csv('vcfQcMetr.tsv', sep = '\t')
-  annot = read.csv('sampleAnnotations.tsv', sep = '\t')
-  samplepc = read.csv('samplePCs.tsv', sep = '\t')
-  refpc = read.csv('refPCs.tsv', sep ='\t')
-  stratify = c('ANCESTRY', 'SeqTech')
-  sds = sampleDataset(bamQcMetrInput = bamQcMetr, vcfQcMetrInput = vcfQcMetr,
-                      annot = annot, primaryID = 'SampleID')
-  return(sds)
-}
diff --git a/R/flagSamples.R b/R/flagSamples.R
index 374f158..ffa0f0f 100644
--- a/R/flagSamples.R
+++ b/R/flagSamples.R
@@ -66,12 +66,12 @@ flagSamples.sampleDataset <- function(sds, cutoffs, zscore = NULL){
     stop("None of metrics in the cutoff table is in the SampleDataset,
          please double check your input data.")
   }
-  if (!is.null(zscore) && !('zscore' %in% attributes(sds)$names)) {
+  if (!is.null(zscore) && !('zscore' %in% names(attributes(sds)))) {
     stop("Input dataset has no z-score, please calculate z-score before
          applying z-score filters.")
   }
 
-  sds$df$flaggedReason = ''
+  sds$df$flaggedReason <- ''
   for (i in 1:dim(cutoffs)[1]) {
     # to do: re-write with sapply
     sds$df = flagSamples(sds$df, cutoffs$qcMetrics[i], cutoffs$value[i],
@@ -84,6 +84,6 @@ flagSamples.sampleDataset <- function(sds, cutoffs, zscore = NULL){
       sds$df = flagSamples(sds$df, i, -zscore, greater = F)
     }
   }
-  sds['flaggedReason'] = 'flaggedReason'
+  sds['flaggedReason'] <- 'flaggedReason'
   return(sds)
 }
diff --git a/R/plot.R b/R/plot.R
index 94dafc3..dd6b843 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -19,27 +19,24 @@ sampleQcPlot.default <- function(
   main = 'QC', geom = c('scatter', 'violin', 'hist')
 ) {
   geom <- match.arg(geom)
-  if (!is.null(annot)) {
-    if (geom == 'scatter') {
-      # scatter plot stratified by sample
-      plt = scatter(data = data, x = 'index', y = qcMetric, strat=annot,
-                     xlab = 'samples', legend = legend, main = main,
-                     outliers = outliers, primaryID = sds$primaryID)
-    }
-    if (geom == 'violin') {
-      data[[annot]] = factor(.toSameLength(data[[annot]]))
-      plt <- (ggplot2::ggplot(data, ggplot2::aes_string(annot, qcMetric,
+  if (!is.null(annot) && geom == 'scatter') {
+    # scatter plot stratified by sample
+    plt <- scatter(data = data, x = 'index', y = qcMetric, strat=annot,
+                  xlab = 'samples', legend = legend, main = main,
+                  outliers = outliers, primaryID = primaryID)
+  }else if(!is.null(annot) && geom == 'violin'){
+    data[[annot]] = factor(.toSameLength(data[[annot]]))
+    plt <- (ggplot2::ggplot(data, ggplot2::aes_string(annot, qcMetric,
                                                       color = annot))
-              + ggplot2::geom_violin()
-              + ggplot2::geom_jitter(height = 0, width = 0.3))
-    }
-  } else {
-    if (geom == 'hist')  {
-      plt = ggplot2::qplot(data[[qcMetric]], geom = 'histogram', bins = 100,
-                           xlab = qcMetric)
-    }
+            + ggplot2::geom_violin()
+            + ggplot2::geom_jitter(height = 0, width = 0.3))
+  }else if(geom == 'hist'){
+    plt <- ggplot2::qplot(data[[qcMetric]], geom = 'histogram', bins = 100,
+                         xlab = qcMetric)
+  }
+  if(!is.null(plt)){
+    plt <- plt + ggplot2::ggtitle(main)
   }
-  plt = plt + ggplot2::ggtitle(main)
   return(plt)
 }
 
@@ -65,17 +62,17 @@ sampleQcPlot.sampleDataset <- function(
   ncols = 5, show = FALSE, sort = TRUE
 ) {
   if(length(qcMetrics) == 1) { ncols = 1 }
-  if (!is.null(outliers)) {
-    if(!all(outliers %in% sds$df[[sds$primaryID]])) {
-      stop("Not all outliers are in the Sample Dataset. Please double check.")
-    }
+  if (!is.null(outliers) && length(setdiff(outliers, sds$df[[sds$primaryID]])) > 0) {
+    stop("Not all outliers are in the Sample Dataset. Please double check.")
   }
-  geom <- match.arg(geom)
-  position <- match.arg(position)
   if(is.null(qcMetrics)) {
-    qcMetrics = sds$qcMetrics
+    qcMetrics <- sds$qcMetrics
   }
-  if (sort) sds = sort(sds, by = annot)
+  if (sort) sds <- sort(sds, by = annot)
+
+  geom <- match.arg(geom)
+  position <- match.arg(position)
+
   plots = sapply(
     qcMetrics,
     function(x) sampleQcPlot(
@@ -102,16 +99,16 @@ outlierPlots <- function(...) UseMethod('outlierPlots')
 outlierPlots.default <- function(tab, qcMetrics, strat, main, outliers,
                          primaryID=NULL, type='violin'){
   plots <- list()
-  plots[[1]] = scatter(tab, x = 'sample', y = qcMetrics, strat = strat,
+  plots[[1]] <- scatter(tab, x = 'sample', y = qcMetrics, strat = strat,
                         xlab = 'sample', outliers = outliers, main = main,
                         legend = T, primaryID = primaryID)
   if (type == 'density') {
-    plots[[2]] = (ggplot2::qplot(tab[, qcMetrics], color = tab[, strat],
+    plots[[2]] <- (ggplot2::qplot(tab[, qcMetrics], color = tab[, strat],
                                 geom = "density", main= main, xlab = qcMetrics)
                   + geom_vline(linetype="dashed",
                                xintercept = tab[outlier.index, qcMetrics]))
   } else {
-    plots[[2]] = (ggplot2::ggplot(data = tab,
+    plots[[2]] <- (ggplot2::ggplot(data = tab,
                                   ggplot2::aes_string(x=strat, y=qcMetrics,
                                                       color = strat))
                   + ggplot2::geom_violin()
@@ -119,12 +116,12 @@ outlierPlots.default <- function(tab, qcMetrics, strat, main, outliers,
                   + ggplot2::ggtitle(main)
                   + ggplot2::theme(axis.text.x = ggplot2::element_blank()))
     if(!is.null(outliers)) {
-      plots[[2]] = (plots[[2]]
+      plots[[2]] <- (plots[[2]]
                     + ggplot2::geom_point(data=tab[tab$sampleId %in% outliers,],
                                           colour = 'black', size = 3))
     }
   }
-  multiplot(plotlist = plots, ncols = 2)
+  return(multiplot(plotlist = plots, ncols = 2))
 }
 
 #' Produce outlier plots for a sampleDataset
diff --git a/R/runExplorer.R b/R/runExplorer.R
index b696812..7ac45d9 100644
--- a/R/runExplorer.R
+++ b/R/runExplorer.R
@@ -3,7 +3,7 @@
 runQcExplorer <- function () {
   appDir <- system.file("apps", "qcExplorer", package = "samplyzer")
   if (appDir == "") {
-    stop("Could not find example directory. Try re-installing `samplyzer`.",
+    stop("Could not find qcExplorer directory. Try re-installing `samplyzer`.",
          call. = FALSE)
   }
   shiny::runApp(appDir, display.mode = "normal")
diff --git a/R/sampleDataset.R b/R/sampleDataset.R
index 425ab8b..aa55521 100644
--- a/R/sampleDataset.R
+++ b/R/sampleDataset.R
@@ -1,8 +1,8 @@
 .loadInput <- function (input) {
   if (is.data.frame(input)) {
-    df = input
+    df <- input
   } else if (is.character(input)) {
-    df = read.csv(input, sep='\t')
+    df <- read.csv(input, sep='\t')
   } else {
     stop("Input is not a tsv file or a data.frame")
   }
@@ -10,12 +10,13 @@
 }
 
 .mergeMetr <- function(df, metrDf, primaryID) {
-  df = merge(df, metrDf, by = primaryID)
+  df <- merge(df, metrDf, by = primaryID)
   return(df)
 }
 
 .getNames <- function(df, primaryID) {
-  df = names(df)[-which(names(df) == primaryID)]
+  return(names(df)[-which(names(df) == primaryID)])
+  #df = names(df)[-which(names(df) == primaryID)]
 }
 
 #' Sample dataset object
@@ -53,37 +54,38 @@ sampleDataset <- function(
 
   # import annotation
   if (is.null(df)) {
-    df = .loadInput(annotInput)
-    annotNames = .getNames(df, primaryID)
+    df <- .loadInput(annotInput)
+    annotNames <- .getNames(df, primaryID)
 
     if (!all(stratify %in% annotNames)) {
       stop("Not all stratified factors are in the annotation file,
            please double check your input stratify option.")
     }
-    bamQcMetrName = NULL
-    vcfQcMetrName = NULL
+
+    bamQcMetrName <- NULL
+    vcfQcMetrName <- NULL
     if (!is.null(bamQcInput)) {
-      bamQcMetr = .loadInput(bamQcInput)
-      df = .mergeMetr(df, bamQcMetr, primaryID)
-      bamQcMetrName = .getNames(bamQcMetr, primaryID)
+      bamQcMetr <- .loadInput(bamQcInput)
+      df <- .mergeMetr(df, bamQcMetr, primaryID)
+      bamQcMetrName <- .getNames(bamQcMetr, primaryID)
     }
     if (!is.null(vcfQcInput)) {
-      vcfQcMetr = .loadInput(vcfQcInput)
-      df = .mergeMetr(df, vcfQcMetr, primaryID)
-      vcfQcMetrName = .getNames(vcfQcMetr, primaryID)
+      vcfQcMetr <- .loadInput(vcfQcInput)
+      df <- .mergeMetr(df, vcfQcMetr, primaryID)
+      vcfQcMetrName <- .getNames(vcfQcMetr, primaryID)
     }
   } else {
     # directly create SDS from an already defined data.frame
-    bamQcMetrName = bamQcInput
-    vcfQcMetrName = vcfQcInput
-    annot = annotInput
+    bamQcMetrName <- bamQcInput
+    vcfQcMetrName <- vcfQcInput
+    annot <- annotInput
   }
 
   sds = list(df = df, qcMetrics = c(bamQcMetrName, vcfQcMetrName),
              bamQcMetr = bamQcMetrName, vcfQcMetr = vcfQcMetrName,
              annot = annotNames, primaryID = primaryID)
 
-  sds$df$index = 1:length(df[[primaryID]])
+  sds$df$index <- 1:length(df[[primaryID]])
   class(sds) <- 'sampleDataset'
   return(sds)
 }