diff --git a/R/gset-fisher.r b/R/gset-fisher.r
index f54fed96..dd0b35e1 100644
--- a/R/gset-fisher.r
+++ b/R/gset-fisher.r
@@ -191,7 +191,8 @@ gset.fisher <- function(genes, genesets, background = NULL,
     d1 <- d + 1 * (d == 0) ## hack to avoid crash...
     b1 <- b + 1 * (b == 0) ## hack to avoid crash...
     pv1 <- try(
-      corpora::fisher.pval(a[ii], (a + b1)[ii], c[ii], (c + d1)[ii], alternative = "greater"), silent = TRUE
+      corpora::fisher.pval(a[ii], (a + b1)[ii], c[ii], (c + d1)[ii], alternative = "greater"),
+      silent = TRUE
     )
     if (class(pv1) != "try-error") {
       pv[ii] <- pv1
diff --git a/R/ngs-fit.r b/R/ngs-fit.r
index a5859a8e..68d43aac 100644
--- a/R/ngs-fit.r
+++ b/R/ngs-fit.r
@@ -59,25 +59,25 @@ ngs.fitContrastsWithAllMethods <- function(counts, X = NULL, samples, design, co
                                            conform.output = TRUE, do.filter = TRUE,
                                            remove.batch = TRUE,
                                            methods = c(
-                                               "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
-                                               "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
+                                             "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
+                                             "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
                                            ),
                                            correct.AveExpr = TRUE, custom = NULL, custom.name = NULL) {
-    ## --------------------------------------------------------------
-    ## Run all tests on raw counts
-    ## --------------------------------------------------------------
-
-    ## Do not test features with full missingness.
-    ## Put them back in the TopTable
-    ## counts0 <- counts
-    ## nas <- apply(counts, 1, function(x) sum(is.na(x)))
-    ## Ex <- names(nas)[which(nas == ncol(counts))]
-    ## keep <- names(nas)[which(nas != ncol(counts))]
-    ## if (length(Ex) > 0) counts <- counts[keep, ]
-  counts <- counts[which(rowMeans(is.na(counts)) < 1),  ]
-  
+  ## --------------------------------------------------------------
+  ## Run all tests on raw counts
+  ## --------------------------------------------------------------
+
+  ## Do not test features with full missingness.
+  ## Put them back in the TopTable
+  ## counts0 <- counts
+  ## nas <- apply(counts, 1, function(x) sum(is.na(x)))
+  ## Ex <- names(nas)[which(nas == ncol(counts))]
+  ## keep <- names(nas)[which(nas != ncol(counts))]
+  ## if (length(Ex) > 0) counts <- counts[keep, ]
+  counts <- counts[which(rowMeans(is.na(counts)) < 1), ]
+
   if (!is.null(X)) {
-    X <- X[which(rowMeans(is.na(X)) < 1),  ]
+    X <- X[which(rowMeans(is.na(X)) < 1), ]
     ## X0 <- X
     ## nas <- apply(X, 1, function(x) sum(is.na(x)))
     ## Ex <- names(nas)[which(nas == ncol(X))]
@@ -85,392 +85,392 @@ ngs.fitContrastsWithAllMethods <- function(counts, X = NULL, samples, design, co
     ## if (length(Ex) > 0) X <- X[keep, ]
   }
 
-    if (methods[1] == "*") {
-        methods <- c(
-            "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
-            "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
-        )
-    }
-    methods <- intersect(methods, c(
-                                      "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
-                                      "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
-                                  ))
-
-    ## If degenerate set design to NULL
-    if (!is.null(design) && ncol(design) >= ncol(X)) {
-        ## "no-replicate" design!!!
-        cat("WARNING: degenerate design. setting design to NULL\n")
-        contr.matrix <- design %*% contr.matrix
-        design <- NULL
-    }
+  if (methods[1] == "*") {
+    methods <- c(
+      "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
+      "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
+    )
+  }
+  methods <- intersect(methods, c(
+    "ttest", "ttest.welch", "voom.limma", "trend.limma", "notrend.limma",
+    "deseq2.wald", "deseq2.lrt", "edger.qlf", "edger.lrt"
+  ))
+
+  ## If degenerate set design to NULL
+  if (!is.null(design) && ncol(design) >= ncol(X)) {
+    ## "no-replicate" design!!!
+    cat("WARNING: degenerate design. setting design to NULL\n")
+    contr.matrix <- design %*% contr.matrix
+    design <- NULL
+  }
 
-    ## ------------------------------------------------------------------
-    ## define transformation methods: log2CPM for counts
-    ## ------------------------------------------------------------------
-    if (is.null(X)) {
-        message("[ngs.fitContrastsWithAllMethods] prior CPM counts =", prior.cpm)
-        message("[ngs.fitContrastsWithAllMethods] CPM scale =", cpm.scale)
-        X <- log2(t(t(counts) / Matrix::colSums(counts)) * cpm.scale + prior.cpm) ## log2CPM
-        X <- limma::normalizeQuantiles(X) ## in linear space
-    } else {
-        message("[ngs.fitContrastsWithAllMethods] using input log-expression matrix X as is")
-    }
+  ## ------------------------------------------------------------------
+  ## define transformation methods: log2CPM for counts
+  ## ------------------------------------------------------------------
+  if (is.null(X)) {
+    message("[ngs.fitContrastsWithAllMethods] prior CPM counts =", prior.cpm)
+    message("[ngs.fitContrastsWithAllMethods] CPM scale =", cpm.scale)
+    X <- log2(t(t(counts) / Matrix::colSums(counts)) * cpm.scale + prior.cpm) ## log2CPM
+    X <- limma::normalizeQuantiles(X) ## in linear space
+  } else {
+    message("[ngs.fitContrastsWithAllMethods] using input log-expression matrix X as is")
+  }
 
-    ## ------------------------------------------------------------------
-    ## get main grouping variable for modeling
-    ## ------------------------------------------------------------------
-    group <- NULL
-    ## if(all(rownames(contr.matrix) %in% samples$group)) {
-    ## }
-    if (!is.null(design)) {
-        group <- colnames(design)[max.col(design)]
-        if (nrow(design) == ncol(design) &&
-            all(rownames(design) == colnames(design))) {
-            group <- NULL
-        }
+  ## ------------------------------------------------------------------
+  ## get main grouping variable for modeling
+  ## ------------------------------------------------------------------
+  group <- NULL
+  ## if(all(rownames(contr.matrix) %in% samples$group)) {
+  ## }
+  if (!is.null(design)) {
+    group <- colnames(design)[max.col(design)]
+    if (nrow(design) == ncol(design) &&
+      all(rownames(design) == colnames(design))) {
+      group <- NULL
     }
+  }
 
-    timings <- list()
-    outputs <- list()
-
-    ## Skip tests that do not tolerate missing values. Inform the user.
-    nmissing <- sum(is.na(X))  
-    
-    ## ---------------- t-Test methods -------------------
-    if ("ttest" %in% methods) {
-        message("[ngs.fitContrastsWithAllMethods] fitting using Student's t-test")
-        timings[["ttest"]] <- system.time(
-            outputs[["ttest"]] <- ngs.fitContrastsWithTTEST(
-                X, contr.matrix, design,
-                method = "equalvar",
-                conform.output = conform.output
-            )
-        )
-    }
+  timings <- list()
+  outputs <- list()
 
-    if ("ttest.rank" %in% methods) {
-        message("[ngs.fitContrastsWithAllMethods] fitting using t-test on ranks")
-        rX <- scale(apply(X, 2, rank, na.last = "keep"))
-        timings[["ttest.rank"]] <- system.time(
-            outputs[["ttest.rank"]] <- ngs.fitContrastsWithTTEST(
-                rX, contr.matrix, design,
-                method = "equalvar",
-                conform.output = conform.output
-            )
-        )
-    }
+  ## Skip tests that do not tolerate missing values. Inform the user.
+  nmissing <- sum(is.na(X))
 
-    if ("ttest.welch" %in% methods) {
-        message("[ngs.fitContrastsWithAllMethods] fitting using Welch t-test")
-        timings[["ttest.welch"]] <- system.time(
-            outputs[["ttest.welch"]] <- ngs.fitContrastsWithTTEST(
-                X, contr.matrix, design,
-                method = "welch",
-                conform.output = conform.output
-            )
-        )
-    }
+  ## ---------------- t-Test methods -------------------
+  if ("ttest" %in% methods) {
+    message("[ngs.fitContrastsWithAllMethods] fitting using Student's t-test")
+    timings[["ttest"]] <- system.time(
+      outputs[["ttest"]] <- ngs.fitContrastsWithTTEST(
+        X, contr.matrix, design,
+        method = "equalvar",
+        conform.output = conform.output
+      )
+    )
+  }
 
-    ## ---------------- LIMMA methods -------------------
-    if ("trend.limma" %in% methods) {
-        message("[ngs.fitContrastsWithAllMethods] fitting using LIMMA trend")
-        tt <- system.time(
-            outputs[["trend.limma"]] <- ngs.fitContrastsWithLIMMA(
-                X, contr.matrix, design,
-                method = "limma", trend = TRUE,
-                prune.samples = prune.samples,
-                conform.output = conform.output, plot = FALSE
-            )
-        )
-        timings[["trend.limma"]] <- round(as.numeric(tt), digits = 4)
-    }
-    if (TRUE && "notrend.limma" %in% methods) {
-        message("[ngs.fitContrastsWithAllMethods] fitting using LIMMA no-trend")
-        timings[["notrend.limma"]] <- system.time(
-            outputs[["notrend.limma"]] <- ngs.fitContrastsWithLIMMA(
-                X, contr.matrix, design,
-                method = "limma", trend = FALSE,
-                prune.samples = prune.samples,
-                conform.output = conform.output, plot = FALSE
-            )
-        )
-    }
-    if ("voom.limma" %in% methods) {
-        if (nmissing>0) {
-            message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
+  if ("ttest.rank" %in% methods) {
+    message("[ngs.fitContrastsWithAllMethods] fitting using t-test on ranks")
+    rX <- scale(apply(X, 2, rank, na.last = "keep"))
+    timings[["ttest.rank"]] <- system.time(
+      outputs[["ttest.rank"]] <- ngs.fitContrastsWithTTEST(
+        rX, contr.matrix, design,
+        method = "equalvar",
+        conform.output = conform.output
+      )
+    )
+  }
+
+  if ("ttest.welch" %in% methods) {
+    message("[ngs.fitContrastsWithAllMethods] fitting using Welch t-test")
+    timings[["ttest.welch"]] <- system.time(
+      outputs[["ttest.welch"]] <- ngs.fitContrastsWithTTEST(
+        X, contr.matrix, design,
+        method = "welch",
+        conform.output = conform.output
+      )
+    )
+  }
+
+  ## ---------------- LIMMA methods -------------------
+  if ("trend.limma" %in% methods) {
+    message("[ngs.fitContrastsWithAllMethods] fitting using LIMMA trend")
+    tt <- system.time(
+      outputs[["trend.limma"]] <- ngs.fitContrastsWithLIMMA(
+        X, contr.matrix, design,
+        method = "limma", trend = TRUE,
+        prune.samples = prune.samples,
+        conform.output = conform.output, plot = FALSE
+      )
+    )
+    timings[["trend.limma"]] <- round(as.numeric(tt), digits = 4)
+  }
+  if (TRUE && "notrend.limma" %in% methods) {
+    message("[ngs.fitContrastsWithAllMethods] fitting using LIMMA no-trend")
+    timings[["notrend.limma"]] <- system.time(
+      outputs[["notrend.limma"]] <- ngs.fitContrastsWithLIMMA(
+        X, contr.matrix, design,
+        method = "limma", trend = FALSE,
+        prune.samples = prune.samples,
+        conform.output = conform.output, plot = FALSE
+      )
+    )
+  }
+  if ("voom.limma" %in% methods) {
+    if (nmissing > 0) {
+      message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
                                                   Cannot perform DGE testing with voom/LIMMA.
                                                   Skipping. If you wish to use voom/LIMMA,
                                                   please impute the data earlier in the platform.")
-        } else {
-            message("[ngs.fitContrastsWithAllMethods] fitting using voom/LIMMA ")
-            timings[["voom.limma"]] <- system.time(
-                outputs[["voom.limma"]] <- ngs.fitContrastsWithLIMMA(
-                    X, contr.matrix, design,
-                    method = "voom",
-                    prune.samples = prune.samples,
-                    conform.output = conform.output, plot = FALSE
-                )
-            )
-        }
+    } else {
+      message("[ngs.fitContrastsWithAllMethods] fitting using voom/LIMMA ")
+      timings[["voom.limma"]] <- system.time(
+        outputs[["voom.limma"]] <- ngs.fitContrastsWithLIMMA(
+          X, contr.matrix, design,
+          method = "voom",
+          prune.samples = prune.samples,
+          conform.output = conform.output, plot = FALSE
+        )
+      )
     }
+  }
 
-    ## ---------------- EdgeR methods -------------------
-    if ("edger.qlf" %in% methods) {
-        if (nmissing>0) {
-            message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
+  ## ---------------- EdgeR methods -------------------
+  if ("edger.qlf" %in% methods) {
+    if (nmissing > 0) {
+      message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
                                                   Cannot perform DGE testing with edgeR using QL F-test.
                                                   Skipping. If you wish to use edgeR with QL F-test,
                                                   please impute the data earlier in the platform.")
-        } else {
-            message("[ngs.fitContrastsWithAllMethods] fitting edgeR using QL F-test ")
-            timings[["edger.qlf"]] <- system.time(
-                outputs[["edger.qlf"]] <- ngs.fitContrastsWithEDGER(
-                    counts, group, contr.matrix, design,
-                    method = "qlf", X = X,
-                    prune.samples = prune.samples,
-                    conform.output = conform.output, plot = FALSE
-                )
-            )
-        }
+    } else {
+      message("[ngs.fitContrastsWithAllMethods] fitting edgeR using QL F-test ")
+      timings[["edger.qlf"]] <- system.time(
+        outputs[["edger.qlf"]] <- ngs.fitContrastsWithEDGER(
+          counts, group, contr.matrix, design,
+          method = "qlf", X = X,
+          prune.samples = prune.samples,
+          conform.output = conform.output, plot = FALSE
+        )
+      )
     }
+  }
 
-    if ("edger.lrt" %in% methods) {
-        if (nmissing>0) {
-            message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
+  if ("edger.lrt" %in% methods) {
+    if (nmissing > 0) {
+      message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
                                                   Cannot perform DGE testing with edgeR using LRT.
                                                   Skipping. If you wish to use edgeR with LRT,
                                                   please impute the data earlier in the platform.")
-        } else {
-            message("[ngs.fitContrastsWithAllMethods] fitting edgeR using LRT")
-            timings[["edger.lrt"]] <- system.time(
-                outputs[["edger.lrt"]] <- ngs.fitContrastsWithEDGER(
-                    counts, group, contr.matrix, design,
-                    method = "lrt", X = X,
-                    prune.samples = prune.samples,
-                    conform.output = conform.output, plot = FALSE
-                )
-            )
-        }
+    } else {
+      message("[ngs.fitContrastsWithAllMethods] fitting edgeR using LRT")
+      timings[["edger.lrt"]] <- system.time(
+        outputs[["edger.lrt"]] <- ngs.fitContrastsWithEDGER(
+          counts, group, contr.matrix, design,
+          method = "lrt", X = X,
+          prune.samples = prune.samples,
+          conform.output = conform.output, plot = FALSE
+        )
+      )
     }
+  }
 
-    ## ---------------- DESEQ2 methods -------------------
-    if ("deseq2.wald" %in% methods) {
-        if (nmissing>0) {
-            message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
+  ## ---------------- DESEQ2 methods -------------------
+  if ("deseq2.wald" %in% methods) {
+    if (nmissing > 0) {
+      message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
                                                   Cannot perform DGE testing with DESeq2 Wald test.
                                                   Skipping. If you wish to use DESeq2 Wald test,
                                                   please impute the data earlier in the platform.")
-        } else {
-            message("[ngs.fitContrastsWithAllMethods] fitting using DESeq2 (Wald test)")
-            timings[["deseq2.wald"]] <- system.time(
-                outputs[["deseq2.wald"]] <- ngs.fitConstrastsWithDESEQ2(
-                    counts, group, contr.matrix, design,
-                    X = X, genes = genes,
-                    test = "Wald", prune.samples = prune.samples,
-                    conform.output = conform.output
-                )
-            )
-        }
+    } else {
+      message("[ngs.fitContrastsWithAllMethods] fitting using DESeq2 (Wald test)")
+      timings[["deseq2.wald"]] <- system.time(
+        outputs[["deseq2.wald"]] <- ngs.fitConstrastsWithDESEQ2(
+          counts, group, contr.matrix, design,
+          X = X, genes = genes,
+          test = "Wald", prune.samples = prune.samples,
+          conform.output = conform.output
+        )
+      )
     }
+  }
 
-    if ("deseq2.lrt" %in% methods) {
-        if (nmissing>0) {
-            message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
+  if ("deseq2.lrt" %in% methods) {
+    if (nmissing > 0) {
+      message("[ngs.fitContrastsWithAllMethods] Detected missing values (NAs) in the data.
                                                   Cannot perform DGE testing with DESeq2 LRT.
                                                   Skipping. If you wish to use DESeq2 LRT,
                                                   please impute the data earlier in the platform.")
-        } else {
-            message("[ngs.fitContrastsWithAllMethods] fitting using DESeq2 (LRT test)")
-            timings[["deseq2.lrt"]] <- system.time(
-                outputs[["deseq2.lrt"]] <- ngs.fitConstrastsWithDESEQ2(
-                    counts, group, contr.matrix, design,
-                    X = X, genes = genes,
-                    test = "LRT", prune.samples = prune.samples,
-                    conform.output = conform.output
-                )
-            )
-        }
+    } else {
+      message("[ngs.fitContrastsWithAllMethods] fitting using DESeq2 (LRT test)")
+      timings[["deseq2.lrt"]] <- system.time(
+        outputs[["deseq2.lrt"]] <- ngs.fitConstrastsWithDESEQ2(
+          counts, group, contr.matrix, design,
+          X = X, genes = genes,
+          test = "LRT", prune.samples = prune.samples,
+          conform.output = conform.output
+        )
+      )
     }
+  }
 
-    if (!is.null(custom)) {
-        message("[ngs.fitContrastsWithAllMethods] adding custom results table")
-        if (is.null(custom.name)) custom.name <- "custom"
-        if (!all(c("tables", "expr") %in% names(custom))) {
-            stop("custom must have 'tables' and 'expr'")
-        }
-        need.tests <- names(outputs[[1]]$tables)
-
-        if (!all(need.tests %in% names(custom$tables))) {
-            stop("custom must include tables: ", paste(need.tests, collapse = " "))
-        }
-        need.cols <- c("external_gene_name", "AveExpr", "adj.P.Val", "P.Value", "logFC")
-        if (!all(need.cols %in% names(custom$tables[[1]]))) {
-            stop("custom tables must include columns: ", paste(need.cols, collapse = " "))
-        }
-        outputs[[custom.name]] <- custom
+  if (!is.null(custom)) {
+    message("[ngs.fitContrastsWithAllMethods] adding custom results table")
+    if (is.null(custom.name)) custom.name <- "custom"
+    if (!all(c("tables", "expr") %in% names(custom))) {
+      stop("custom must have 'tables' and 'expr'")
     }
+    need.tests <- names(outputs[[1]]$tables)
 
-    ## ----------------------------------------------------------------------
-    ## "corrections" ...
-    ## ----------------------------------------------------------------------
-    if (correct.AveExpr) {
-        message("[ngs.fitContrastsWithAllMethods] correcting AveExpr values...")
-        message("[ngs.fitContrastsWithAllMethods] dim.X: ", dim(X)[1], ",", dim(X)[2])
-
-        ## Some methods like edgeR and Deseq2 compute some weird
-        ## normalized expression matrix. We need to "correct" for
-        ## those.
-
-        exp.matrix <- contr.matrix
-        if (!is.null(design)) exp.matrix <- (design %*% contr.matrix)
-        samplesX <- lapply(apply(exp.matrix != 0, 2, which), function(i) rownames(exp.matrix)[i])
-        samples1 <- lapply(apply(exp.matrix > 0, 2, which), function(i) rownames(exp.matrix)[i])
-        samples0 <- lapply(apply(exp.matrix < 0, 2, which), function(i) rownames(exp.matrix)[i])
-
-        avgX <- sapply(samplesX, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
-        avg.1 <- sapply(samples1, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
-        avg.0 <- sapply(samples0, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
-
-        dim(avgX)
-        i <- j <- 1
-        for (i in 1:length(outputs)) {
-            for (j in 1:length(outputs[[i]]$tables)) {
-                outputs[[i]]$tables[[j]]$AveExpr <- avgX[, j]
-                outputs[[i]]$tables[[j]]$AveExpr1 <- avg.1[, j]
-                outputs[[i]]$tables[[j]]$AveExpr0 <- avg.0[, j]
-            }
-        }
+    if (!all(need.tests %in% names(custom$tables))) {
+      stop("custom must include tables: ", paste(need.tests, collapse = " "))
     }
+    need.cols <- c("external_gene_name", "AveExpr", "adj.P.Val", "P.Value", "logFC")
+    if (!all(need.cols %in% names(custom$tables[[1]]))) {
+      stop("custom tables must include columns: ", paste(need.cols, collapse = " "))
+    }
+    outputs[[custom.name]] <- custom
+  }
 
-    ## ----------------------------------------------------------------------
-    ## add some statistics
-    ## ----------------------------------------------------------------------
-    message("[ngs.fitContrastsWithAllMethods] calculating statistics...")
-    i <- 1
+  ## ----------------------------------------------------------------------
+  ## "corrections" ...
+  ## ----------------------------------------------------------------------
+  if (correct.AveExpr) {
+    message("[ngs.fitContrastsWithAllMethods] correcting AveExpr values...")
+    message("[ngs.fitContrastsWithAllMethods] dim.X: ", dim(X)[1], ",", dim(X)[2])
+
+    ## Some methods like edgeR and Deseq2 compute some weird
+    ## normalized expression matrix. We need to "correct" for
+    ## those.
+
+    exp.matrix <- contr.matrix
+    if (!is.null(design)) exp.matrix <- (design %*% contr.matrix)
+    samplesX <- lapply(apply(exp.matrix != 0, 2, which), function(i) rownames(exp.matrix)[i])
+    samples1 <- lapply(apply(exp.matrix > 0, 2, which), function(i) rownames(exp.matrix)[i])
+    samples0 <- lapply(apply(exp.matrix < 0, 2, which), function(i) rownames(exp.matrix)[i])
+
+    avgX <- sapply(samplesX, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
+    avg.1 <- sapply(samples1, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
+    avg.0 <- sapply(samples0, function(s) rowMeans(X[, s, drop = FALSE], na.rm = TRUE))
+
+    dim(avgX)
+    i <- j <- 1
     for (i in 1:length(outputs)) {
-        res <- outputs[[i]]
-        M <- sapply(res$tables, function(x) x[, "AveExpr"]) ## !!!! edgeR and Deseq2 are weird!!!
-        M0 <- sapply(res$tables, function(x) x[, "AveExpr0"])
-        M1 <- sapply(res$tables, function(x) x[, "AveExpr1"])
-        Q <- sapply(res$tables, function(x) x[, "adj.P.Val"])
-        P <- sapply(res$tables, function(x) x[, "P.Value"])
-        logFC <- sapply(res$tables, function(x) x[, "logFC"])
-        colnames(M) <- colnames(logFC) <- colnames(P) <- colnames(Q) <- colnames(contr.matrix)
-        rownames(M) <- rownames(logFC) <- rownames(P) <- rownames(Q) <- rownames(res$tables[[1]])
-        rownames(M0) <- rownames(M1) <- rownames(res$tables[[1]])
-
-        ## count significant terms
-        qvalues <- c(1e-16, 10**seq(-8, -2, 2), 0.05, 0.1, 0.2, 0.5, 1)
-        lfc <- 1
-        sig.both <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (abs(logFC) > lfc), na.rm = TRUE))
-        sig.up <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (logFC > lfc), na.rm = TRUE))
-        sig.down <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (logFC < -lfc), na.rm = TRUE))
-        sig.notsig <- sapply(qvalues, function(q) Matrix::colSums(Q > q | (abs(logFC) < lfc), na.rm = TRUE))
-        if (NCOL(Q) == 1) {
-            sig.both <- matrix(sig.both, nrow = 1)
-            sig.up <- matrix(sig.up, nrow = 1)
-            sig.down <- matrix(sig.down, nrow = 1)
-            sig.notsig <- matrix(sig.notsig, nrow = 1)
-            rownames(sig.both) <- rownames(sig.up) <- colnames(Q)
-            rownames(sig.down) <- rownames(sig.notsig) <- colnames(Q)
-        }
-        colnames(sig.both) <- colnames(sig.notsig) <- qvalues
-        colnames(sig.up) <- colnames(sig.down) <- qvalues
-
-        res$sig.counts <- list(both = sig.both, up = sig.up, down = sig.down, notsig = sig.notsig)
-
-        ## need this? takes space!!!
-        res$p.value <- P
-        res$q.value <- Q
-        res$logFC <- logFC
-        res$aveExpr <- M
-        res$aveExpr0 <- M0
-        res$aveExpr1 <- M1
+      for (j in 1:length(outputs[[i]]$tables)) {
+        outputs[[i]]$tables[[j]]$AveExpr <- avgX[, j]
+        outputs[[i]]$tables[[j]]$AveExpr1 <- avg.1[, j]
+        outputs[[i]]$tables[[j]]$AveExpr0 <- avg.0[, j]
+      }
+    }
+  }
 
-        outputs[[i]] <- res
+  ## ----------------------------------------------------------------------
+  ## add some statistics
+  ## ----------------------------------------------------------------------
+  message("[ngs.fitContrastsWithAllMethods] calculating statistics...")
+  i <- 1
+  for (i in 1:length(outputs)) {
+    res <- outputs[[i]]
+    M <- sapply(res$tables, function(x) x[, "AveExpr"]) ## !!!! edgeR and Deseq2 are weird!!!
+    M0 <- sapply(res$tables, function(x) x[, "AveExpr0"])
+    M1 <- sapply(res$tables, function(x) x[, "AveExpr1"])
+    Q <- sapply(res$tables, function(x) x[, "adj.P.Val"])
+    P <- sapply(res$tables, function(x) x[, "P.Value"])
+    logFC <- sapply(res$tables, function(x) x[, "logFC"])
+    colnames(M) <- colnames(logFC) <- colnames(P) <- colnames(Q) <- colnames(contr.matrix)
+    rownames(M) <- rownames(logFC) <- rownames(P) <- rownames(Q) <- rownames(res$tables[[1]])
+    rownames(M0) <- rownames(M1) <- rownames(res$tables[[1]])
+
+    ## count significant terms
+    qvalues <- c(1e-16, 10**seq(-8, -2, 2), 0.05, 0.1, 0.2, 0.5, 1)
+    lfc <- 1
+    sig.both <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (abs(logFC) > lfc), na.rm = TRUE))
+    sig.up <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (logFC > lfc), na.rm = TRUE))
+    sig.down <- sapply(qvalues, function(q) Matrix::colSums((Q <= q) * (logFC < -lfc), na.rm = TRUE))
+    sig.notsig <- sapply(qvalues, function(q) Matrix::colSums(Q > q | (abs(logFC) < lfc), na.rm = TRUE))
+    if (NCOL(Q) == 1) {
+      sig.both <- matrix(sig.both, nrow = 1)
+      sig.up <- matrix(sig.up, nrow = 1)
+      sig.down <- matrix(sig.down, nrow = 1)
+      sig.notsig <- matrix(sig.notsig, nrow = 1)
+      rownames(sig.both) <- rownames(sig.up) <- colnames(Q)
+      rownames(sig.down) <- rownames(sig.notsig) <- colnames(Q)
     }
+    colnames(sig.both) <- colnames(sig.notsig) <- qvalues
+    colnames(sig.up) <- colnames(sig.down) <- qvalues
 
-    ## --------------------------------------------------------------
-    ## Reshape matrices by comparison
-    ## --------------------------------------------------------------
-    message("[ngs.fitContrastsWithAllMethods] reshape matrices...")
-
-    tests <- colnames(outputs[[1]]$p.value)
-    ntest <- length(tests)
-    P <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$p.value[, i]))
-    Q <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$q.value[, i]))
-    logFC <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$logFC[, i]))
-    aveExpr <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr[, i]))
-    aveExpr0 <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr0[, i]))
-    aveExpr1 <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr1[, i]))
-
-    sig.up <- lapply(1:ntest, function(i) {
-        do.call(rbind, lapply(outputs, function(x) x$sig.counts[["up"]][i, ]))
-    })
-    sig.down <- lapply(1:ntest, function(i) {
-        do.call(rbind, lapply(outputs, function(x) x$sig.counts[["down"]][i, ]))
-    })
-    sig.notsig <- lapply(1:ntest, function(i) {
-        do.call(rbind, lapply(outputs, function(x) x$sig.counts[["notsig"]][i, ]))
-    })
-    sig.both <- lapply(1:ntest, function(i) {
-        do.call(rbind, lapply(outputs, function(x) x$sig.counts[["both"]][i, ]))
-    })
-    names(P) <- names(Q) <- names(logFC) <- names(aveExpr) <- tests
-    names(sig.up) <- names(sig.down) <- names(sig.both) <- names(sig.notsig) <- tests
-    sig.counts <- list(up = sig.up, down = sig.down, both = sig.both, notsig = sig.notsig)
+    res$sig.counts <- list(both = sig.both, up = sig.up, down = sig.down, notsig = sig.notsig)
 
-    ## --------------------------------------------------
-    ## meta analysis, aggregate p-values
-    ## --------------------------------------------------
-    message("[ngs.fitContrastsWithAllMethods] aggregating p-values...")
+    ## need this? takes space!!!
+    res$p.value <- P
+    res$q.value <- Q
+    res$logFC <- logFC
+    res$aveExpr <- M
+    res$aveExpr0 <- M0
+    res$aveExpr1 <- M1
 
-    all.meta <- list()
-    i <- 1
-    for (i in 1:ntest) {
-        pv <- P[[i]]
-        qv <- Q[[i]]
-        fc <- logFC[[i]]
-        mx <- aveExpr[[i]]
-        mx0 <- aveExpr0[[i]]
-        mx1 <- aveExpr1[[i]]
-
-        ## avoid strange values
-        fc[is.infinite(fc) | is.nan(fc)] <- NA
-        pv <- pmax(pv, 1e-99)
-        pv[is.na(pv)] <- 1
-        qv[is.na(qv)] <- 1
-
-
-        ## !!!!!!!!!!!!!!!!!!!!!!!! NEED RETHINK !!!!!!!!!!!!!!!!!!!!!!!!
-        meta.p <- apply(pv, 1, function(p) exp(mean(log(p), na.rm = TRUE))) ## geometric mean
-        meta.q <- stats::p.adjust(meta.p, method = "BH")
-        meta.fx <- rowMeans(fc, na.rm = TRUE)
-        meta.avg <- rowMeans(mx, na.rm = TRUE)
-        meta.avg0 <- rowMeans(mx0, na.rm = TRUE)
-        meta.avg1 <- rowMeans(mx1, na.rm = TRUE)
-
-        meta <- data.frame(fx = meta.fx, p = meta.p, q = meta.q)
-        avg <- data.frame(avg.0 = meta.avg0, avg.1 = meta.avg1)
-        rownames(meta) <- rownames(logFC[[i]])
-        rownames(avg) <- rownames(logFC[[i]])
-        rownames(fc) <- NULL ## saves memory
-        rownames(pv) <- NULL
-        rownames(qv) <- NULL
-        all.meta[[i]] <- data.frame(meta = meta, avg, fc = I(fc), p = I(pv), q = I(qv))
-        if (!is.null(genes)) all.meta[[i]] <- cbind(genes, all.meta[[i]])
-    }
-    names(all.meta) <- tests
+    outputs[[i]] <- res
+  }
 
-    timings0 <- do.call(rbind, timings)
-    colnames(timings0) <- names(timings[[1]])
-    colnames(timings0) <- c("user.self", "sys.self", "elapsed", "user.child", "sys.child")
+  ## --------------------------------------------------------------
+  ## Reshape matrices by comparison
+  ## --------------------------------------------------------------
+  message("[ngs.fitContrastsWithAllMethods] reshape matrices...")
+
+  tests <- colnames(outputs[[1]]$p.value)
+  ntest <- length(tests)
+  P <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$p.value[, i]))
+  Q <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$q.value[, i]))
+  logFC <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$logFC[, i]))
+  aveExpr <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr[, i]))
+  aveExpr0 <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr0[, i]))
+  aveExpr1 <- lapply(1:ntest, function(i) sapply(outputs, function(x) x$aveExpr1[, i]))
+
+  sig.up <- lapply(1:ntest, function(i) {
+    do.call(rbind, lapply(outputs, function(x) x$sig.counts[["up"]][i, ]))
+  })
+  sig.down <- lapply(1:ntest, function(i) {
+    do.call(rbind, lapply(outputs, function(x) x$sig.counts[["down"]][i, ]))
+  })
+  sig.notsig <- lapply(1:ntest, function(i) {
+    do.call(rbind, lapply(outputs, function(x) x$sig.counts[["notsig"]][i, ]))
+  })
+  sig.both <- lapply(1:ntest, function(i) {
+    do.call(rbind, lapply(outputs, function(x) x$sig.counts[["both"]][i, ]))
+  })
+  names(P) <- names(Q) <- names(logFC) <- names(aveExpr) <- tests
+  names(sig.up) <- names(sig.down) <- names(sig.both) <- names(sig.notsig) <- tests
+  sig.counts <- list(up = sig.up, down = sig.down, both = sig.both, notsig = sig.notsig)
+
+  ## --------------------------------------------------
+  ## meta analysis, aggregate p-values
+  ## --------------------------------------------------
+  message("[ngs.fitContrastsWithAllMethods] aggregating p-values...")
+
+  all.meta <- list()
+  i <- 1
+  for (i in 1:ntest) {
+    pv <- P[[i]]
+    qv <- Q[[i]]
+    fc <- logFC[[i]]
+    mx <- aveExpr[[i]]
+    mx0 <- aveExpr0[[i]]
+    mx1 <- aveExpr1[[i]]
+
+    ## avoid strange values
+    fc[is.infinite(fc) | is.nan(fc)] <- NA
+    pv <- pmax(pv, 1e-99)
+    pv[is.na(pv)] <- 1
+    qv[is.na(qv)] <- 1
+
+
+    ## !!!!!!!!!!!!!!!!!!!!!!!! NEED RETHINK !!!!!!!!!!!!!!!!!!!!!!!!
+    meta.p <- apply(pv, 1, function(p) exp(mean(log(p), na.rm = TRUE))) ## geometric mean
+    meta.q <- stats::p.adjust(meta.p, method = "BH")
+    meta.fx <- rowMeans(fc, na.rm = TRUE)
+    meta.avg <- rowMeans(mx, na.rm = TRUE)
+    meta.avg0 <- rowMeans(mx0, na.rm = TRUE)
+    meta.avg1 <- rowMeans(mx1, na.rm = TRUE)
+
+    meta <- data.frame(fx = meta.fx, p = meta.p, q = meta.q)
+    avg <- data.frame(avg.0 = meta.avg0, avg.1 = meta.avg1)
+    rownames(meta) <- rownames(logFC[[i]])
+    rownames(avg) <- rownames(logFC[[i]])
+    rownames(fc) <- NULL ## saves memory
+    rownames(pv) <- NULL
+    rownames(qv) <- NULL
+    all.meta[[i]] <- data.frame(meta = meta, avg, fc = I(fc), p = I(pv), q = I(qv))
+    if (!is.null(genes)) all.meta[[i]] <- cbind(genes, all.meta[[i]])
+  }
+  names(all.meta) <- tests
 
-    res <- list(
-        outputs = outputs, meta = all.meta, sig.counts = sig.counts,
-        timings = timings0, X = X
-    )
-    return(res)
+  timings0 <- do.call(rbind, timings)
+  colnames(timings0) <- names(timings[[1]])
+  colnames(timings0) <- c("user.self", "sys.self", "elapsed", "user.child", "sys.child")
+
+  res <- list(
+    outputs = outputs, meta = all.meta, sig.counts = sig.counts,
+    timings = timings0, X = X
+  )
+  return(res)
 }
 
 ## --------------------------------------------------------------------------------------------
diff --git a/R/pgx-correct.R b/R/pgx-correct.R
index 3e2a6bde..436545c0 100644
--- a/R/pgx-correct.R
+++ b/R/pgx-correct.R
@@ -1398,7 +1398,9 @@ bc.evaluateResults <- function(xlist, pheno, lfc = 0.2, q = 0.2, pos = NULL,
       if (clust == "tsne") {
         nb <- max(0.33, min(30, round(ncol(xlist[[1]]) / 5)))
         ## CLUSTFUN <- function(x) uwot::tumap(scale(t(x), scale = FALSE), n_neighbors = nb)
-        if (ncol(xlist[[1]]) <= 6) { nb <- 0.5 }
+        if (ncol(xlist[[1]]) <= 6) {
+          nb <- 0.5
+        }
         CLUSTFUN <- function(x) {
           Rtsne::Rtsne(scale(t(x)),
             check_duplicates = FALSE,
@@ -1754,7 +1756,9 @@ compare_batchcorrection_methods <- function(X, samples, pheno, contrasts,
   if (clust.method == "tsne" && nmissing == 0) {
     message("Computing t-SNE clustering...")
     nb <- max(0.33, round(min(30, dim(X) / 5)))
-    if (ncol(X) <= 6) { nb <- 0.5 }
+    if (ncol(X) <= 6) {
+      nb <- 0.5
+    }
     pos <- lapply(xlist, function(x) {
       Rtsne::Rtsne(t2(x), perplexity = nb, check_duplicates = FALSE)$Y
     })
@@ -1964,8 +1968,8 @@ svaCorrect <- function(X, y) {
   #  pp <- paste0(model.par, collapse = "+")
   #  lm.expr <- paste0("lm(t(X) ~ ", pp, ", data=pheno)")
   X.r <- t(stats::resid(lm(t(X) ~ y)))
-  ##n.sv <- isva::EstDimRMT(X.r, FALSE)$dim + 1
-  n.sv <- isva::EstDimRMT(X.r, FALSE)$dim 
+  ## n.sv <- isva::EstDimRMT(X.r, FALSE)$dim + 1
+  n.sv <- isva::EstDimRMT(X.r, FALSE)$dim
   n.sv
   ## top 1000 genes only (faster)
   X1 <- Matrix::head(X[order(-matrixStats::rowSds(X, na.rm = TRUE)), ], 1000)
@@ -2075,8 +2079,8 @@ ruvCorrect <- function(X, y, k = NULL, type = c("III", "g"), controls = 0.10) {
   if (any(is.na(X))) {
     stop("[ruvCorrect] cannot handle missing values in X.")
   }
-  sdx <- matrixStats::rowSds(X,na.rm=TRUE)
-  
+  sdx <- matrixStats::rowSds(X, na.rm = TRUE)
+
   ## F-test using limma just variables
   if (!is.null(y) && length(controls) == 1 && is.numeric(controls[1])) {
     ii <- which(!duplicated(rownames(X)) & sdx > 0)
diff --git a/R/pgx-correlation.R b/R/pgx-correlation.R
index d9a8019f..7fe319f5 100644
--- a/R/pgx-correlation.R
+++ b/R/pgx-correlation.R
@@ -400,7 +400,6 @@ pgx.computePartialCorrelationMatrix <- function(tX, method = PCOR.METHODS, fast
 #'
 #' @export
 pgx.testPhenoCorrelation <- function(df, plot = TRUE, cex = 1, compute.pv = TRUE) {
-
   cl <- sapply(df, class)
   nlev <- apply(df, 2, function(x) length(unique(x[!is.na(x)])))
   cvar <- which(cl %in% c("numeric", "integer") & nlev >= 2)
@@ -451,7 +450,7 @@ pgx.testPhenoCorrelation <- function(df, plot = TRUE, cex = 1, compute.pv = TRUE
       rownames(fisher.P) <- colnames(dd)
       colnames(fisher.P) <- colnames(dd)
     }
-    
+
     ## discrete vs continuous -> ANOVA or Kruskal-Wallace
     kruskal.P <- NULL
     if (ncol(dc) > 0) {
@@ -466,7 +465,7 @@ pgx.testPhenoCorrelation <- function(df, plot = TRUE, cex = 1, compute.pv = TRUE
       rownames(kruskal.P) <- colnames(dd)
       colnames(kruskal.P) <- colnames(dc)
     }
-    
+
     ## continuous vs continuous -> correlation test
     cor.P <- NULL
     if (ncol(dc) > 1) {
diff --git a/R/pgx-pathways.R b/R/pgx-pathways.R
index 7269058b..76b9b9fd 100644
--- a/R/pgx-pathways.R
+++ b/R/pgx-pathways.R
@@ -1,25 +1,25 @@
 #' Download and color nodes in PathBank SVGs
 #'
-#' Downloads a PathBank SVG file, modifies its XML structure to color pathway nodes 
+#' Downloads a PathBank SVG file, modifies its XML structure to color pathway nodes
 #' according to input values, and returns the path to the modified file.
 #'
 #' @param pb A character string specifying the PathBank pathway ID.
-#' @param val A named numeric vector where names correspond to pathway element IDs, 
-#'            and values represent the level of regulation (e.g., fold-change or intensity). 
-#'            Positive values indicate up-regulation, while negative values indicate down-regulation. 
+#' @param val A named numeric vector where names correspond to pathway element IDs,
+#'            and values represent the level of regulation (e.g., fold-change or intensity).
+#'            Positive values indicate up-regulation, while negative values indicate down-regulation.
 #'
-#' @return A character string representing the path to the modified SVG file, 
+#' @return A character string representing the path to the modified SVG file,
 #'         or `NULL` if the download fails.
 #'
 #' @details
-#' The function downloads the SVG representation of a PathBank pathway, modifies its 
+#' The function downloads the SVG representation of a PathBank pathway, modifies its
 #' XML content to color pathway nodes based on the provided `val` parameter:
 #' - **Positive values**: Nodes are colored red, with the intensity increasing as the value grows.
 #' - **Negative values**: Nodes are colored blue, with the intensity increasing as the absolute value grows.
 #'
 #' Opacity is determined by the square root of the absolute value, scaled and capped at a maximum of 1.
 #'
-#' Nodes with IDs matching the names in `val` are colored accordingly. Nodes without matching IDs 
+#' Nodes with IDs matching the names in `val` are colored accordingly. Nodes without matching IDs
 #' or `val` values are left unchanged.
 #'
 #' @examples
@@ -55,10 +55,17 @@ pathbankview <- function(pb, val) {
   writeLines(lines, destfile)
 
   # Load the SVG
-  doc <- tryCatch({
-    read_xml(destfile)
-  }, error = function(w) {NULL})
-  if (is.null(doc)) return(NULL)
+  doc <- tryCatch(
+    {
+      read_xml(destfile)
+    },
+    error = function(w) {
+      NULL
+    }
+  )
+  if (is.null(doc)) {
+    return(NULL)
+  }
 
   label_nodes <- xml_find_all(doc, ".//text")
   labels <- xml2::xml_text(label_nodes)
@@ -103,26 +110,26 @@ pathbankview <- function(pb, val) {
 
 #' Download and color nodes in WikiPathways SVGs
 #'
-#' Downloads a WikiPathways SVG file, modifies its XML structure to color pathway nodes 
+#' Downloads a WikiPathways SVG file, modifies its XML structure to color pathway nodes
 #' according to regulation values, and returns the path to the modified file.
 #'
 #' @param wp A character string specifying the WikiPathways pathway ID.
-#' @param val A named numeric vector where names correspond to pathway element IDs, 
-#'            and values represent the level of regulation (e.g., fold-change or intensity). 
-#'            Positive values indicate up-regulation, while negative values indicate down-regulation. 
+#' @param val A named numeric vector where names correspond to pathway element IDs,
+#'            and values represent the level of regulation (e.g., fold-change or intensity).
+#'            Positive values indicate up-regulation, while negative values indicate down-regulation.
 #'
-#' @return A character string representing the path to the modified SVG file, 
+#' @return A character string representing the path to the modified SVG file,
 #'         or `NULL` if the download fails.
 #'
 #' @details
-#' The function downloads the SVG representation of a WikiPathways pathway, modifies its 
+#' The function downloads the SVG representation of a WikiPathways pathway, modifies its
 #' XML content to color pathway nodes based on the provided `val` parameter:
 #' - **Positive values**: Nodes are colored red, with the intensity increasing as the value grows.
 #' - **Negative values**: Nodes are colored blue, with the intensity increasing as the absolute value grows.
 #'
 #' Opacity is determined by the square root of the absolute value, scaled and capped at a maximum of 1.
 #'
-#' Nodes with IDs matching the names in `val` are colored accordingly. Nodes without matching IDs 
+#' Nodes with IDs matching the names in `val` are colored accordingly. Nodes without matching IDs
 #' or `val` values are left unchanged.
 #' @examples
 #' \dontrun{
@@ -178,10 +185,17 @@ wikipathview <- function(wp, val) {
   writeLines(lines, destfile)
 
   # Load the SVG
-  doc <- tryCatch({
-    read_xml(destfile)
-  }, error = function(w) {NULL})
-  if (is.null(doc)) return(NULL)
+  doc <- tryCatch(
+    {
+      read_xml(destfile)
+    },
+    error = function(w) {
+      NULL
+    }
+  )
+  if (is.null(doc)) {
+    return(NULL)
+  }
 
   # Find all 'text' elements
   label_nodes <- xml_find_all(doc, ".//text")
diff --git a/R/pgx-pheno.R b/R/pgx-pheno.R
index 79b32469..32dbd9be 100644
--- a/R/pgx-pheno.R
+++ b/R/pgx-pheno.R
@@ -396,7 +396,7 @@ expandPhenoMatrix <- function(M, drop.ref = TRUE, keep.numeric = FALSE, check =
   y.isnum <- (y.class %in% c("numeric", "integer"))
   kk <- which(y.isnum | (!y.isnum & nlevel > 1 & nratio < 0.66))
   if (length(kk) == 0) {
-    kk <- which(y.isnum | (!y.isnum & nlevel > 1 ))
+    kk <- which(y.isnum | (!y.isnum & nlevel > 1))
   }
   if (length(kk) == 0) {
     return(NULL)
diff --git a/tests/testthat/test-pgx-pathways.R b/tests/testthat/test-pgx-pathways.R
index f079ba25..4b75f2a8 100644
--- a/tests/testthat/test-pgx-pathways.R
+++ b/tests/testthat/test-pgx-pathways.R
@@ -1,12 +1,11 @@
 test_that("pathbankview works correctly", {
-
   # Test 1: Valid PathBank ID and coloring node
   test_that("Valid PathBank ID with coloring", {
     pb <- "SMP0080852"
     val <- c("PW_C000064" = 2.0)
     result <- pathbankview(pb, val)
-    expect_type(result, "character")  # Expect a character string (path to SVG)
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string (path to SVG)
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 2: Valid PathBank ID with no matching nodes
@@ -14,8 +13,8 @@ test_that("pathbankview works correctly", {
     pb <- "SMP0080852"
     val <- c("nonexistent_node" = 1.5)
     result <- pathbankview(pb, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 3: Invalid PathBank ID
@@ -23,7 +22,7 @@ test_that("pathbankview works correctly", {
     pb <- "INVALID_ID"
     val <- c("PW_C000064" = 2.0)
     result <- pathbankview(pb, val)
-    expect_null(result)  # Expect NULL when the ID is invalid
+    expect_null(result) # Expect NULL when the ID is invalid
   })
 
   # Test 4: Valid PathBank ID with NULL val
@@ -31,8 +30,8 @@ test_that("pathbankview works correctly", {
     pb <- "SMP0080852"
     val <- NULL
     result <- pathbankview(pb, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 5: Valid PathBank ID with mixed values in val
@@ -40,8 +39,8 @@ test_that("pathbankview works correctly", {
     pb <- "SMP0080852"
     val <- c("PW_C000064" = 1.5, "PW_C000065" = -2.0)
     result <- pathbankview(pb, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 6: Empty val vector
@@ -49,20 +48,19 @@ test_that("pathbankview works correctly", {
     pb <- "SMP0080852"
     val <- c()
     result <- pathbankview(pb, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 })
 
 test_that("wikipath (modern) works correctly", {
-
   # Test 1: Valid WikiPathway ID and coloring node
   test_that("Valid WikiPathway ID with coloring", {
     wp <- "WP5405"
     val <- c("CRYL1" = 2.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string (path to SVG)
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string (path to SVG)
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 2: Valid WikiPathway ID with no matching nodes
@@ -70,8 +68,8 @@ test_that("wikipath (modern) works correctly", {
     wp <- "WP5405"
     val <- c("node" = 2.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 3: Invalid WikiPathway ID
@@ -79,7 +77,7 @@ test_that("wikipath (modern) works correctly", {
     wp <- "invalid"
     val <- c("CRYL1" = 2.0)
     result <- wikipathview(wp, val)
-    expect_null(result)  # Expect NULL when the ID is invalid
+    expect_null(result) # Expect NULL when the ID is invalid
   })
 
   # Test 4: Valid WikiPathway ID with NULL val
@@ -87,8 +85,8 @@ test_that("wikipath (modern) works correctly", {
     wp <- "WP5405"
     val <- NULL
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 5: Valid WikiPathway ID with mixed values in val
@@ -96,8 +94,8 @@ test_that("wikipath (modern) works correctly", {
     wp <- "WP5405"
     val <- c("CRYL1" = 2.0, "CRYL2" = 1.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 6: Empty val vector
@@ -105,20 +103,19 @@ test_that("wikipath (modern) works correctly", {
     wp <- "WP5405"
     val <- c()
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 })
 
 test_that("wikipath (classic) works correctly", {
-
   # Test 1: Valid WikiPathway ID and coloring node
   test_that("Valid WikiPathway ID with coloring", {
     wp <- "WP4876"
     val <- c("TRAF3" = 2.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string (path to SVG)
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string (path to SVG)
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 2: Valid WikiPathway ID with no matching nodes
@@ -126,8 +123,8 @@ test_that("wikipath (classic) works correctly", {
     wp <- "WP4876"
     val <- c("node" = 2.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 3: Invalid WikiPathway ID
@@ -135,7 +132,7 @@ test_that("wikipath (classic) works correctly", {
     wp <- "invalid"
     val <- c("TRAF3" = 2.0)
     result <- wikipathview(wp, val)
-    expect_null(result)  # Expect NULL when the ID is invalid
+    expect_null(result) # Expect NULL when the ID is invalid
   })
 
   # Test 4: Valid WikiPathway ID with NULL val
@@ -143,8 +140,8 @@ test_that("wikipath (classic) works correctly", {
     wp <- "WP4876"
     val <- NULL
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 5: Valid WikiPathway ID with mixed values in val
@@ -152,8 +149,8 @@ test_that("wikipath (classic) works correctly", {
     wp <- "WP4876"
     val <- c("TRAF3" = 2.0, "TRAF4" = 1.0)
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 
   # Test 6: Empty val vector
@@ -161,7 +158,7 @@ test_that("wikipath (classic) works correctly", {
     wp <- "WP4876"
     val <- c()
     result <- wikipathview(wp, val)
-    expect_type(result, "character")  # Expect a character string
-    expect_true(file.exists(result))  # Ensure the returned path exists
+    expect_type(result, "character") # Expect a character string
+    expect_true(file.exists(result)) # Ensure the returned path exists
   })
 })