@@ -1,6 +1,6 @@

 Package: GRaNIE

 Title: GRaNIE: Reconstruction cell type specific gene regulatory networks including enhancers using chromatin accessibility and RNA-seq data

-Version: 1.1.14

+Version: 1.1.15

 Encoding: UTF-8

 Authors@R: c(person("Christian", "Arnold", email =

         "chrarnold@web.de", role = c("cre","aut")),

@@ -3,7 +3,7 @@

 #' Genetic variants associated with diseases often affect non-coding regions, thus likely having a regulatory role. To understand the effects of genetic variants in these regulatory regions, identifying genes that are modulated by specific regulatory elements (REs) is crucial. The effect of gene regulatory elements, such as enhancers, is often cell-type specific, likely because the combinations of transcription factors (TFs) that are regulating a given enhancer have celltype specific activity. This TF activity can be quantified with existing tools such as \code{diffTF} and captures differences in binding of a TF in open chromatin regions. Collectively, this forms a gene regulatory network (\code{eGRN}) with cell-type and data-specific TF-RE and RE-gene links. Here, we reconstruct such a \code{eGRN} using bulk RNAseq and open chromatin (e.g., using ATACseq or ChIPseq for open chromatin marks) and optionally TF activity data. Our network contains different types of links, connecting TFs to regulatory elements, the latter of which is connected to genes in the vicinity or within the same chromatin domain (TAD). We use a statistical framework to assign empirical FDRs and weights to all links using a permutation-based approach.

 #' 

 #' @section Package functions:

-#' See the Vignettes for a workflow example and more generally \url{https://grp-zaugg.embl-community.io/GRaNIE/articles/} for all project-related information.

+#' See the Vignettes for a workflow example and more generally \url{https://grp-zaugg.embl-community.io/GRaNIE} for all project-related information.

 #' 

 #' @section GRN object:

 #' The \code{GRaNIE} package works with \code{GRN} objects. See \code{\linkS4class{GRN}} for details.

@@ -39,8 +39,14 @@ initializeGRN <- function(objectMetadata = list(),

   .checkAndLoadPackagesGenomeAssembly(genomeAssembly)

   # Create the folder first if not yet existing

+  checkmate::assertCharacter(outputFolder, min.chars = 1, len = 1)

   if (!dir.exists(outputFolder)) {

-    dir.create(outputFolder)

+    res = dir.create(outputFolder)

+    if (!res) {

+        message = paste0("Could not create the output directory ", outputFolder, ". Check the path /and/or access rights.")

+        .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

+    }

+    checkmate::assertDirectoryExists(outputFolder, access = "w")

   }

   # Create an absolute path out of the given outputFolder now that it exists

@@ -109,6 +115,9 @@ initializeGRN <- function(objectMetadata = list(),

   futile.logger::flog.info(paste0("Empty GRN object created successfully. Type the object name (e.g., GRN) to retrieve summary information about it at any time."))

+  futile.logger::flog.info(paste0(" Default output folder: ", GRN@config$directories$outputRoot))

+  futile.logger::flog.info(paste0(" Genome assembly: ", genomeAssembly))

+  

   .printExecutionTime(start, prefix = "")

   GRN

@@ -187,11 +196,17 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     # Normalize ID column names

     if (idColumn_peaks != "peakID") {

       counts_peaks = dplyr::rename(counts_peaks, peakID = !!(idColumn_peaks))

+      idColumn_peaks ="peakID"

     }

     if (idColumn_RNA != "ENSEMBL") {

       counts_rna = dplyr::rename(counts_rna, ENSEMBL = !!(idColumn_RNA))

+      idColumn_RNA = "ENSEMBL"

     }

-    

+      

+    # Check existence of correct ID column now

+    checkmate::assertSubset(idColumn_peaks, colnames(counts_peaks))

+    checkmate::assertSubset(idColumn_RNA, colnames(counts_rna))

+

     # Check ID columns for missing values and remove

     rna_missing_ID =  which(is.na(counts_rna$ENSEMBL))

     if (length(rna_missing_ID) > 0) {

@@ -257,8 +272,11 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     GRN@data$peaks$counts_norm = data.l[["peaks"]] %>% dplyr::mutate(isFiltered = FALSE)

     countsRNA.norm.df = data.l[["RNA"]]  %>% dplyr::mutate(isFiltered = FALSE)

+    futile.logger::flog.info(paste0( "Final dimensions of data:"))

+    futile.logger::flog.info(paste0( " RNA  : ", nrow(countsRNA.norm.df)  , " x ", ncol(countsRNA.norm.df)   - 2, " (rows x columns)"))

+    futile.logger::flog.info(paste0( " peaks: ", nrow(countsPeaks.norm.df), " x ", ncol(countsPeaks.norm.df) - 1, " (rows x columns)"))

     # Create permutations for RNA

-    futile.logger::flog.info(paste0( "Produce ", .getMaxPermutation(GRN), " permutations of RNA-counts"))

+    futile.logger::flog.info(paste0( "Generate ", .getMaxPermutation(GRN), " permutations of RNA-counts"))

     GRN@data$RNA$counts_norm.l = .shuffleColumns(countsRNA.norm.df, .getMaxPermutation(GRN), returnUnshuffled = TRUE, returnAsList = TRUE)

     if (!is.null(sampleMetadata)) {

@@ -341,7 +359,7 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     overlappingPeaks = which(GenomicRanges::countOverlaps(consensus.gr ,consensus.gr) >1)

-    if (length(overlappingPeaks) > 0){

+    if (length(overlappingPeaks) > 0) {

       ids = (consensus.gr[overlappingPeaks] %>% as.data.frame())$peakID

@@ -574,11 +592,12 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 }

 .populateGeneAnnotation <- function (GRN) {

-  

-  futile.logger::flog.info(paste0(" Calculate statistics for each gene (mean and CV)"))

-  

+

+

   countsRNA.clean  = getCounts(GRN, type = "rna", norm = TRUE, permuted = FALSE)

+  futile.logger::flog.info(paste0(" Calculate statistics for each of the ", nrow(countsRNA.clean), " genes that were provided with the RNA-seq data (mean and CV)"))

+  

   countsRNA.m = as.matrix(dplyr::select(countsRNA.clean, -.data$ENSEMBL))

   rowMeans_rna = rowMeans(countsRNA.m)

@@ -591,7 +610,7 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

                                 gene.mean = rowMeans_rna, 

                                 gene.median = rowMedians_rna, 

                                 gene.CV = CV_rna) %>%

-    dplyr::full_join(genomeAnnotation.df, by = c("gene.ENSEMBL")) %>%

+    dplyr::left_join(genomeAnnotation.df, by = c("gene.ENSEMBL")) %>%

     dplyr::mutate(gene.type = forcats::fct_explicit_na(.data$gene.type, na_level = "unknown/missing"))

   GRN@annotation$genes = metadata_rna

@@ -704,6 +723,8 @@ filterData <- function (GRN,

   checkmate::assertNumber(maxNormalizedMean_peaks, lower = minNormalizedMean_peaks , null.ok = TRUE)

   checkmate::assertNumber(maxNormalizedMeanRNA, lower = minNormalizedMeanRNA, null.ok = TRUE)

   checkmate::assertCharacter(chrToKeep_peaks, min.len = 1, any.missing = FALSE, null.ok = TRUE)

+  checkmate::assertSubset(chrToKeep_peaks, GRN@data$peaks$consensusPeaks %>% dplyr::pull(.data$chr) %>% unique() %>% as.character())

+  

   checkmate::assertIntegerish(minSize_peaks, lower = 1, null.ok = TRUE)

   checkmate::assertIntegerish(maxSize_peaks, lower = dplyr::if_else(is.null(minSize_peaks), 1, minSize_peaks), null.ok = TRUE)

   checkmate::assertNumber(minCV_peaks, lower = 0, null.ok = TRUE)

@@ -739,6 +760,11 @@ filterData <- function (GRN,

     futile.logger::flog.info(paste0("Collectively, filter ", nPeaksBefore -length(peaks_toKeep), " out of ", nPeaksBefore, " peaks."))

     futile.logger::flog.info(paste0("Number of remaining peaks: ", length(peaks_toKeep)))

+    if (length(peaks_toKeep) < 1000) {

+        message = paste0("Too few peaks (", length(peaks_toKeep), ") remain after filtering. At least 1000 peaks must remain. Adjust the filtering settings.")

+        .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

+    }

+    

     GRN@data$peaks$consensusPeaks$isFiltered  = ! GRN@data$peaks$consensusPeaks$peakID  %in% peaks_toKeep

     #GRN@data$peaks$counts_raw$isFiltered = ! GRN@data$peaks$counts_raw$peakID  %in% peaks_toKeep

     GRN@data$peaks$counts_norm$isFiltered = ! GRN@data$peaks$counts_norm$peakID  %in% peaks_toKeep

@@ -757,7 +783,12 @@ filterData <- function (GRN,

                                     minMean = minNormalizedMeanRNA, maxMean = maxNormalizedMeanRNA, 

                                     minCV = minCV_genes, maxCV = maxCV_genes) 

-    futile.logger::flog.info(paste0(" Number of rows in total: ", nrow(GRN@data$RNA$counts_norm.l[["0"]])))

+    

+    if (length(genes.CV) < 100) {

+        message = paste0("Too few genes (", length(genes.CV), ") remain after filtering. At least 100 genes must remain. Adjust the filtering settings.")

+        .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

+    }

+    

     for (permutationCur in c(0)) {

       permIndex = as.character(permutationCur)

       rowMeans = rowMeans(dplyr::select(GRN@data$RNA$counts_norm.l[[permIndex]], -.data$ENSEMBL))

@@ -768,8 +799,6 @@ filterData <- function (GRN,

     # Raw counts are left untouched and filtered where needed only

     futile.logger::flog.info(paste0(" Flagged ", nRowsFlagged, " rows because the row mean was smaller than ", minNormalizedMeanRNA))

-    # TODO: Filter genes by CV

-    

     GRN@config$isFiltered = TRUE

   } 

@@ -1093,7 +1122,7 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

   if (is.null(GRN@data$TFs$TF_peak_overlap) | forceRerun) {

-    futile.logger::flog.info(paste0("Overlap peaks and TFBS using ", nCores, " cores. This may take a few minutes..."))

+    futile.logger::flog.info(paste0("Overlap peaks and TFBS using ", nCores, " cores. This may take a while, particularly if the number of samples is large..."))

     genomeAssembly = GRN@config$parameters$genomeAssembly

     seqlengths = .getChrLengths(genomeAssembly)

@@ -1254,17 +1283,25 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

   start = Sys.time()

+  # Set the column name to just ENSEMBL to avoid column mismatch issues

+  HOCOMOCO_mapping$ENSEMBL = HOCOMOCO_mapping$TF.ENSEMBL

+  

   # Filter to only the TFs

   # In addition, the no of TF because multiple TFs can map to the same gene/ ENSEMBL ID

   # Also filter 0 count genes because they otherwise cause errors downstream

   rowSums = rowSums(dplyr::select(matrix1, -.data$ENSEMBL))

   # Keep only Ensembl IDs from TFs we have data from

-  matrix1.norm.TFs.df = dplyr::filter(matrix1, .data$ENSEMBL %in% HOCOMOCO_mapping$ENSEMBL, rowSums != 0)

+  matrix1.norm.TFs.df = dplyr::filter(matrix1, .data$ENSEMBL %in% HOCOMOCO_mapping$TF.ENSEMBL, rowSums != 0)

+  

+  nFiltered1 = dplyr::filter(matrix1, ! .data$ENSEMBL %in% HOCOMOCO_mapping$TF.ENSEMBL) %>% nrow()

+  nFiltered2 = dplyr::filter(matrix1, rowSums == 0) %>% nrow()

   diff = nrow(matrix1) - nrow(matrix1.norm.TFs.df)

   if (diff > 0) {

-    message = paste0(whitespacePrefix, "Retain ", nrow(matrix1.norm.TFs.df), " rows from TF/gene data out of ", nrow(matrix1), " (filter non-TF genes and TF genes with 0 counts throughout and keep only unique ENSEMBL IDs).")

+    message = paste0(whitespacePrefix, "Retain ", nrow(matrix1.norm.TFs.df), " unique genes from TF/gene data out of ", nrow(matrix1), 

+                     " (filter ",  nFiltered1, " non-TF genes and ", nFiltered2, 

+                     " TF genes with 0 counts throughout).")

     futile.logger::flog.info(message)

   }

@@ -1273,7 +1310,7 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

     .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

   }

-  HOCOMOCO_mapping.exp = dplyr::filter(HOCOMOCO_mapping, .data$ENSEMBL %in% matrix1.norm.TFs.df$ENSEMBL)

+  HOCOMOCO_mapping.exp = dplyr::filter(HOCOMOCO_mapping, .data$TF.ENSEMBL %in% matrix1.norm.TFs.df$ENSEMBL)

   futile.logger::flog.info(paste0(whitespacePrefix, "Correlate TF/gene data for ", nrow(matrix1.norm.TFs.df), " unique Ensembl IDs (TFs) and peak counts for ", nrow(matrix_peaks), " peaks."))

   futile.logger::flog.info(paste0(whitespacePrefix, "Note: For subsequent steps, the same gene may be associated with multiple TF, depending on the translation table."))

   # Correlate TF gene counts with peak counts 

@@ -1301,10 +1338,10 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

   # Reorder to make sure the order is the same. Due to the duplication ID issue, the number of columns may increase after the column selection

   # Some columns may be removed here due to zero standard deviation

-  HOCOMOCO_mapping.exp.filt = HOCOMOCO_mapping.exp %>% dplyr::filter(.data$ENSEMBL %in% colnames(sort.cor.m))

+  HOCOMOCO_mapping.exp.filt = HOCOMOCO_mapping.exp %>% dplyr::filter(.data$TF.ENSEMBL %in% colnames(sort.cor.m))

   sort.cor.m = sort.cor.m[,as.character(HOCOMOCO_mapping.exp.filt$ENSEMBL)] 

-  colnames(sort.cor.m) = as.character(HOCOMOCO_mapping.exp.filt$HOCOID)

+  colnames(sort.cor.m) = as.character(HOCOMOCO_mapping.exp.filt$TF.HOCOID)

   .printExecutionTime(start, prefix = whitespacePrefix)

   sort.cor.m

@@ -2007,14 +2044,16 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

                                      corMethod,

                                      whitespacePrefix = "  ")

-    # TODO:Check here for the TFs

+    # TODO:Check here for the TFs BMAL1.0.A

     allTF = intersect(colnames(peak_TF_overlap.df), colnames(peaksCor.m))

+    checkmate::assertIntegerish(length(allTF), lower = 1)

+    

     futile.logger::flog.info(paste0(" Run FDR calculations for ", length(allTF), " TFs for which TFBS predictions and ",

                                     connectionTypeCur, " data for the corresponding gene are available."))

     if (length(allTF) < ncol(peak_TF_overlap.df) | length(allTF) < ncol(peaksCor.m)) {

       TF_missing = setdiff(colnames(peak_TF_overlap.df), colnames(peaksCor.m))

-      if (length(TF_missing) > 0) futile.logger::flog.info(paste0("  Skip the following ", length(TF_missing), " TF due to missing data: ", paste0(TF_missing, collapse = ",")))

+      if (length(TF_missing) > 0) futile.logger::flog.info(paste0("  Skip the following ", length(TF_missing), " TF due to missing dataor because they are marked as filtered: ", paste0(TF_missing, collapse = ",")))

     }

     stopifnot(nrow(peaksCor.m) == nrow(peak_TF_overlap.df))

@@ -160,7 +160,9 @@

     countsRNA.df  = dplyr::select(countsRNA,  tidyselect::all_of(c(idColumn_RNA,  sharedColumns)))

     countsPeaks.df = dplyr::select(countsPeaks, tidyselect::all_of(c(idColumn_peaks, sharedColumns)))

-    futile.logger::flog.info(paste0(" ", length(sharedColumns), " samples (", paste0(sharedColumns, collapse = ","), ") are shared between the peaks and RNA-Seq data"))

+    futile.logger::flog.info(paste0(" ", length(sharedColumns), " samples (", 

+                                    dplyr::if_else(length(sharedColumns) > 100, "names omitted here", paste0(sharedColumns, collapse = ",")), 

+                                    ") are shared between the peaks and RNA-Seq data"))

     notIntersecting_peaks = setdiff(colnames(countsPeaks),c(sharedColumns,idColumn_peaks))

     notIntersecting_RNA  = setdiff(colnames(countsRNA),c(sharedColumns,idColumn_RNA))

@@ -998,6 +998,12 @@ calculateCommunitiesEnrichment <- function(GRN,

       .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

   }

+  if (length(selCommunities) < length(communities)) {

+      missingCommunities = setdiff(communities, selCommunities)

+      message = paste0("Some of the requested communities (", paste0(missingCommunities , collapse = ","), ") were not found. Only the following communities are available: ", paste0(selCommunities, collapse = ","))

+      .checkAndLogWarningsAndErrors(NULL, message, isWarning = TRUE)

+  }

+  

   selCommunities

 }

@@ -326,7 +326,9 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

           plot(g)

         }, error = function(e) {

-          futile.logger::flog.warn(paste0(" Could not plot PCA with variable ", varCur))

+          message = paste0(" Could not plot PCA with variable ", varCur)

+          .checkAndLogWarningsAndErrors(NULL, message, isWarning = TRUE)   

+            

           plot(c(0, 1), c(0, 1), ann = FALSE, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')

           message = paste0(" Could not plot PCA with variable ", varCur)

           text(x = 0.5, y = 0.5, message, cex = 1.6, col = "red")

@@ -493,7 +495,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

 #' @template basenameOutput

 #' @template plotDetails

 #' @param dataType Character vector. One of, or both of, \code{"real"} or \code{"permuted"}. For which data type, real or permuted data, to produce the diagnostic plots?

-#' @param nTFMax \code{NULL} or Integer. Default \code{NULL}. Maximum number of TFs to process. Can be used for testing purposes by setting this to a small number i(.e., 10)

+#' @param nTFMax \code{NULL} or Integer > 0. Default \code{NULL}. Maximum number of TFs to process. Can be used for testing purposes by setting this to a small number i(.e., 10)

 #' @template plotAsPDF

 #' @template pdf_width

 #' @param pdf_height_base  Number. Default 8. Base height of the PDF, in cm, per connection type. The total height is automatically determined based on the number of connection types that are found in the object (e.g., expression or TF activity). For example, when two connection types are found, the base height is multiplied by 2.

@@ -524,7 +526,7 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

   checkmate::assert(checkmate::checkNull(basenameOutput), checkmate::checkCharacter(basenameOutput, len = 1, min.chars = 1, any.missing = FALSE))

   checkmate::assertFlag(plotDetails)

   checkmate::assertSubset(dataType, c("real", "permuted"), empty.ok = FALSE)

-  checkmate::assert(checkmate::checkNull(nTFMax), checkmate::checkIntegerish(nTFMax))

+  checkmate::assert(checkmate::checkNull(nTFMax), checkmate::checkIntegerish(nTFMax, lower = 1))

   checkmate::assertFlag(plotAsPDF)

   checkmate::assertNumeric(pdf_width, lower = 5, upper = 99)

   checkmate::assertNumeric(pdf_height_base, lower = 5, upper = 99)

@@ -917,6 +919,10 @@ plotDiagnosticPlots_peakGene <- function(GRN,

   checkmate::assert(checkmate::checkNull(outputFolder), checkmate::checkCharacter(outputFolder, min.chars = 1))

   checkmate::assert(checkmate::checkNull(basenameOutput), checkmate::checkCharacter(basenameOutput, len = 1, min.chars = 1, any.missing = FALSE))

   checkmate::assertList(gene.types, any.missing = FALSE, min.len = 1, types = "character")

+  for (geneTypesCur in gene.types) {

+      checkmate::assertSubset(geneTypesCur, c("all", unique(as.character(GRN@annotation$genes$gene.type))) %>% stats::na.omit(), empty.ok = FALSE)

+  }  

+

   checkmate::assertFlag(useFiltered)

   checkmate::assertFlag(plotDetails)

   checkmate::assertFlag(plotPerTF)

@@ -1122,6 +1128,9 @@ plotDiagnosticPlots_peakGene <- function(GRN,

         filenameCur = paste0(fileBase, paste0(geneTypesSelected, collapse = "+"), filteredStr, ".pdf")

         .checkOutputFile(filenameCur)

         grDevices::pdf(file = filenameCur, width = pdf_width, height = pdf_height)

+        

+        futile.logger::flog.info(paste0(" Plotting to file ", filenameCur))

+        

       }

       if (plotPerTF) {

@@ -2342,8 +2351,8 @@ plotGeneralEnrichment <- function(GRN, outputFolder = NULL, basenameOutput = NUL

 #' @inheritParams plotCommunitiesEnrichment

 #' @param communities Numeric vector. Default \code{seq_len(10)}. Depending on what was specified in the \code{display} parameter, this parameter would indicate either the rank or the label of the communities to be plotted. i.e. for \code{communities = c(1,4)}, if \code{display = "byRank"} the results for the first and fourth largest communities will be plotted. if \code{display = "byLabel"}, the results for the communities labeled \code{"1"}, and \code{"4"} will be plotted. If set to \code{NULL}, all communities will be plotted

 #' @template forceRerun

-#' @param topnGenes Integer. Default 20. Number of genes to plot, sorted by their rank or label.

-#' @param topnTFs Integer. Default 20. Number of TFs to plot, sorted by their rank or label.

+#' @param topnGenes Integer > 0. Default 20. Number of genes to plot, sorted by their rank or label.

+#' @param topnTFs Integer > 0. Default 20. Number of TFs to plot, sorted by their rank or label.

 #' @return The same \code{\linkS4class{GRN}} object, without modifications.

 #' @seealso \code{\link{plotGeneralGraphStats}}

 #' @seealso \code{\link{calculateCommunitiesStats}}

@@ -2354,7 +2363,7 @@ plotGeneralEnrichment <- function(GRN, outputFolder = NULL, basenameOutput = NUL

 #' GRN = plotCommunitiesStats(GRN, plotAsPDF = FALSE, pages = 1)

 #' @export

 plotCommunitiesStats <- function(GRN, outputFolder = NULL, basenameOutput = NULL, 

-                                 display = "byRank", communities = seq_len(10), 

+                                 display = "byRank", communities = seq_len(5), 

                                  topnGenes = 20, topnTFs = 20, 

                                  plotAsPDF = TRUE, pdf_width = 12, pdf_height = 12, pages = NULL,

                                  forceRerun = FALSE){

@@ -2368,8 +2377,8 @@ plotCommunitiesStats <- function(GRN, outputFolder = NULL, basenameOutput = NULL

   checkmate::assert(checkmate::checkNull(outputFolder), checkmate::checkCharacter(outputFolder, min.chars = 1))

   checkmate::assert(checkmate::checkNull(basenameOutput), checkmate::checkCharacter(basenameOutput, len = 1, min.chars = 1, any.missing = FALSE))

-  checkmate::assertIntegerish(topnGenes)

-  checkmate::assertIntegerish(topnTFs)

+  checkmate::assertIntegerish(topnGenes, lower = 1)

+  checkmate::assertIntegerish(topnTFs, lower = 1)

   checkmate::assertChoice(display , c("byRank", "byLabel"))

   checkmate::assert(checkmate::checkNull(communities), checkmate::checkNumeric(communities, lower = 1, any.missing = FALSE, min.len = 1))

   checkmate::assertFlag(plotAsPDF)

@@ -2542,8 +2551,8 @@ plotCommunitiesStats <- function(GRN, outputFolder = NULL, basenameOutput = NULL

 #' @inheritParams plotGeneralEnrichment

 #' @param display Character. Default \code{"byRank"}. One of: \code{"byRank"}, \code{"byLabel"}. Specify whether the communities will be displayed based on their rank, where the largest community (with most vertices) would have a rank of 1, or by their label. Note that the label is independent of the rank.

 #' @param communities \code{NULL} or numeric vector. Default \code{NULL}. If set to \code{NULL}, the default, all communities enrichments that have been calculated before are plotted. If a numeric vector is specified: Depending on what was specified in the \code{display} parameter, this parameter indicates either the rank or the label of the communities to be plotted. i.e. for \code{communities = c(1,4)}, if \code{display = "byRank"} the results for the first and fourth largest communities are plotted. if \code{display = "byLabel"}, the results for the communities labeled \code{"1"}, and \code{"4"} are plotted. 

-#' @param nSignificant Numeric >0. Default 3. Threshold to filter out an ontology term with less than \code{nSignificant} overlapping genes. 

-#' @param nID Numeric >0. Default 10. For the reduced summary heatmap, number of top terms to select per community / for the general enrichment.

+#' @param nSignificant Numeric > 0. Default 3. Threshold to filter out an ontology term with less than \code{nSignificant} overlapping genes. 

+#' @param nID Numeric > 0. Default 10. For the reduced summary heatmap, number of top terms to select per community / for the general enrichment.

 #' @template maxWidth_nchar_plot

 #' @return  The same \code{\linkS4class{GRN}} object, without modifications.

 #' @seealso \code{\link{plotGeneralEnrichment}}

@@ -2921,8 +2930,8 @@ plotTFEnrichment <- function(GRN, rankType = "degree", n = NULL, TF.names = NULL

   checkmate::assert(checkmate::checkNull(n), checkmate::checkNumeric(n))

   checkmate::assertNumeric(p, lower = 0, upper = 1)

   checkmate::assertNumeric(topn_pvalue, lower = 1)

-  checkmate::assertNumeric(nSignificant, lower = 0)

-  checkmate::assertNumeric(nID, lower = 0)

+  checkmate::assertNumeric(nSignificant, lower = 1)

+  checkmate::assertNumeric(nID, lower = 1)

   checkmate::assertFlag(plotAsPDF)

   checkmate::assertNumeric(pdf_width, lower = 5, upper = 99)

   checkmate::assertNumeric(pdf_height, lower = 5, upper = 99)

@@ -2938,13 +2947,16 @@ plotTFEnrichment <- function(GRN, rankType = "degree", n = NULL, TF.names = NULL

   outputFolder = .checkOutputFolder(GRN, outputFolder)

-  if (rankType == "custom"){

+  if (rankType == "custom") {

     if(is.null(TF.names)){

-      futile.logger::flog.error("To plot the GO enrichment for a custom set of TFs, you must provide the TF names in the 'TF.names' parameter.")

+      message = "To plot the GO enrichment for a custom set of TFs, you must provide the TF names in the 'TF.names' parameter."

+      .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

     }

     wrongTFs = setdiff(TF.names, unique(GRN@connections$all.filtered$`0`$TF.name))

     if (length(wrongTFs)>0){

-      futile.logger::flog.warn(paste0("The TFs ",  paste0(wrongTFs, collapse = " + "), " are not in the filtered GRN. They will be ommited from the results."))

+       message = paste0("The TFs ",  paste0(wrongTFs, collapse = " + "), " are not in the filtered GRN. They will be ommited from the results.")

+       .checkAndLogWarningsAndErrors(NULL, message, isWarning = TRUE)

+

     }

     TFset = setdiff(TF.names, wrongTFs) 

   } else{ #rankType = "degree"

@@ -3467,7 +3479,8 @@ visualizeGRN <- function(GRN, outputFolder = NULL,  basenameOutput = NULL, plotA

     if (nrow(grn.merged) == 0) {

-        futile.logger::flog.warn(paste0("No rows left in the GRN. Creating empty plot."))

+        message = paste0("No rows left in the GRN. Creating empty plot.")

+        .checkAndLogWarningsAndErrors(NULL, message, isWarning = TRUE)

         plot(c(0, 1), c(0, 1), ann = FALSE, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n', main = title)

         message = paste0(title, "\n\nThe GRN has no edges that pass the filter criteria.")

@@ -43,9 +43,9 @@ As they are listed under \code{Suggests}, they may not yet be installed, and the

 \item{display}{Character. Default \code{"byRank"}. One of: \code{"byRank"}, \code{"byLabel"}. Specify whether the communities will be displayed based on their rank, where the largest community (with most vertices) would have a rank of 1, or by their label. Note that the label is independent of the rank.}

-\item{topnGenes}{Integer. Default 20. Number of genes to plot, sorted by their rank or label.}

+\item{topnGenes}{Integer > 0. Default 20. Number of genes to plot, sorted by their rank or label.}

-\item{topnTFs}{Integer. Default 20. Number of TFs to plot, sorted by their rank or label.}

+\item{topnTFs}{Integer > 0. Default 20. Number of TFs to plot, sorted by their rank or label.}

 \item{maxWidth_nchar_plot}{Integer (>=10). Default 50. Maximum number of characters for a term before it is truncated.}

@@ -9,7 +9,7 @@ plotCommunitiesStats(

   outputFolder = NULL,

   basenameOutput = NULL,

   display = "byRank",

-  communities = seq_len(10),

+  communities = seq_len(5),

   topnGenes = 20,

   topnTFs = 20,

   plotAsPDF = TRUE,

@@ -30,9 +30,9 @@ plotCommunitiesStats(

 \item{communities}{Numeric vector. Default \code{seq_len(10)}. Depending on what was specified in the \code{display} parameter, this parameter would indicate either the rank or the label of the communities to be plotted. i.e. for \code{communities = c(1,4)}, if \code{display = "byRank"} the results for the first and fourth largest communities will be plotted. if \code{display = "byLabel"}, the results for the communities labeled \code{"1"}, and \code{"4"} will be plotted. If set to \code{NULL}, all communities will be plotted}

-\item{topnGenes}{Integer. Default 20. Number of genes to plot, sorted by their rank or label.}

+\item{topnGenes}{Integer > 0. Default 20. Number of genes to plot, sorted by their rank or label.}

-\item{topnTFs}{Integer. Default 20. Number of TFs to plot, sorted by their rank or label.}

+\item{topnTFs}{Integer > 0. Default 20. Number of TFs to plot, sorted by their rank or label.}

 \item{plotAsPDF}{\code{TRUE} or \code{FALSE}. Default \code{TRUE}.Should the plots be printed to a PDF file? If set to \code{TRUE}, a PDF file is generated, the name of which depends on the value of \code{basenameOutput}. If set to \code{FALSE}, all plots are printed to the currently active device. Note that most functions print more than one plot, which means you may only see the last plot depending on your active graphics device.}