new file mode 100644

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

+^.*\.Rproj$

+^\.Rproj\.user$

+pkgdown/

+TFBS_selected/

+exampleDataset/

+logo.png

@@ -4,3 +4,14 @@ vignettes/*_files

 vignettes/*_cache

 vignettes/output

 .Rproj.user

+docs

+exampleDataset

+pkgdown

+TFBS_selected

+TFBS_selected.zip

+_pkgdown.yml

+.Rhistory

+.RData

+GRaNIE.Rproj

+fixBugs_RELASE_GRaNIE.sh

+logo.png

@@ -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.10

+Version: 1.1.12

 Encoding: UTF-8

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

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

@@ -27,28 +27,22 @@ Imports:

     RColorBrewer,

     ComplexHeatmap,

     DESeq2,

-    csaw,

     circlize,

-    robust,

     progress,

     utils,

     methods,

     stringr,

     scales,

-    BiocManager,

-    BiocParallel,

     igraph,

     S4Vectors,

     ggplot2,

     rlang,

     Biostrings,

     GenomeInfoDb,

-    IRanges,

     SummarizedExperiment,

     forcats,

     gridExtra,

     limma,

-    purrr,

     tidyselect,

     readr,

     grid,

@@ -60,12 +54,10 @@ Imports:

     magrittr,

     tibble,

     viridis,

-    BiocFileCache,

-    colorspace

+    colorspace,

+    biomaRt

 Depends:

-    R (>= 4.2.0),

-    tidyverse,

-    topGO

+    R (>= 4.2.0)

 Suggests:

     knitr,

     BSgenome.Hsapiens.UCSC.hg19,

@@ -79,13 +71,17 @@ Suggests:

     org.Hs.eg.db,

     org.Mm.eg.db,

     IHW,

-    biomaRt,

     clusterProfiler,

     ReactomePA,

     DOSE,

+    BiocFileCache,

     ChIPseeker,

     testthat (>= 3.0.0),

-    BiocStyle

+    BiocStyle,

+    csaw,

+    BiocParallel,

+    topGO,

+    robust

 VignetteBuilder: knitr

 biocViews: Software, GeneExpression, GeneRegulation, NetworkInference, GeneSetEnrichment, BiomedicalInformatics, Genetics, Transcriptomics, ATACSeq, RNASeq, GraphAndNetwork, Regression, Transcription, ChIPSeq

 License: Artistic-2.0

@@ -24,6 +24,7 @@ export(initializeGRN)

 export(loadExampleObject)

 export(nGenes)

 export(nPeaks)

+export(nTFs)

 export(overlapPeaksAndTFBS)

 export(performAllNetworkAnalyses)

 export(plotCommunitiesEnrichment)

@@ -36,8 +37,6 @@ export(plotPCA_all)

 export(plotTFEnrichment)

 export(plot_stats_connectionSummary)

 export(visualizeGRN)

-import(BiocFileCache)

-import(BiocManager)

 import(GenomicRanges)

 import(checkmate)

 import(ggplot2)

@@ -45,11 +44,11 @@ import(grDevices)

 import(graphics)

 import(patchwork)

 import(tibble)

-import(tidyverse)

-import(topGO)

 import(utils)

 importFrom(GenomicRanges,mcols)

-importFrom(IRanges,width)

+importFrom(Matrix,Matrix)

+importFrom(biomaRt,getBM)

+importFrom(biomaRt,useEnsembl)

 importFrom(circlize,colorRamp2)

 importFrom(grid,gpar)

 importFrom(magrittr,`%>%`)

@@ -1,3 +1,19 @@

+# GRaNIE 1.1.12 (2022-09-13)

+

+## Major changes

+

+- many Documentation and R help updates, Package Details Vignette is online

+- The workflow vignette is now improved: better figure resolution, figure aspect ratios are optimized, and a few other changes

+- the eGRN graph structure as built by `build_eGRN_graph()` in the `GRaNIE` object is now reset whenever the function `filterGRNAndConnectGenes()` is successfully executed to make sure that enrichment functions etc are not using an outdated graph structure. 

+- the landing page of the website has been extended and overhauled

+- removed some dependency packages and moved others into `Suggests` to lower the installation burden of the package. In addition, removed `topGO` from the `Depends` section (now in `Suggests`) and removed `tidyverse` altogether (before in `Depends`). Detailed explanations when and how the packages listed under `Suggests` are needed can now be found in the new [Package Details Vignette](https://grp-zaugg.embl-community.io/GRaNIE/articles/GRaNIE_packageDetails.html#installation).

+- new helper function `installSuggestedPackages()` to install all packages listed under `Suggests`

+

+## Minor changes

+

+- many small changes in the code, better argument checking, and preparing rigorous unit test inclusion

+- internally renaming the (recently changed / renamed) gene type `lncRNA` from `biomaRt` to `lincRNA` to be compatible with older versions of `GRaNIE`

+

 # GRaNIE 1.1.X (2022-05-31)

@@ -202,7 +202,7 @@ setMethod("show",

                 cat(" Communities (TF-gene):\n")

                 df = igraph::vertex.attributes(GRN@graph[["TF_gene"]]$graph)

                 if (!is.null(df) & "community" %in% names(df)) {

-                    communities = df %>% as.data.frame() %>% dplyr::count(community) %>% dplyr::arrange(desc(n))

+                    communities = df %>% as.data.frame() %>% dplyr::count(community) %>% dplyr::arrange(dplyr::desc(.data$n))

                     cat("  Communities, sorted by size (n = Number of nodes): ", paste0(communities$community, " (n=", communities$n, collapse = "), "), ")\n", sep = "")

                 } else {

                     cat("  None found\n")

@@ -241,11 +241,13 @@ setMethod("show",

 #' Get the number of peaks for a \code{\linkS4class{GRN}} object.

 #' 

-#' Return the number of peaks (all or only non-filtered ones) that are defined in the \code{\linkS4class{GRN}} object.

+#' Returns the number of peaks (all or only non-filtered ones) from the provided peak datain the \code{\linkS4class{GRN}} object.

 #'

 #' @template GRN

 #' @param filter TRUE or FALSE. Default TRUE. Should peaks marked as filtered be included in the count?

 #' @return Integer. Number of peaks that are defined in the \code{\linkS4class{GRN}} object, either by excluding (filter = TRUE) or including (filter = FALSE) peaks that are currently marked as \emph{filtered}.

+#' @seealso \code{\link{nTFs}}

+#' @seealso \code{\link{nGenes}}

 #' @examples

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -253,7 +255,6 @@ setMethod("show",

 #' nPeaks(GRN, filter = FALSE)

 #' @export

 #' @aliases peaks

-#' @rdname peaks-methods

 nPeaks <- function(GRN, filter = TRUE) {

   checkmate::assertClass(GRN, "GRN")

@@ -276,11 +277,13 @@ nPeaks <- function(GRN, filter = TRUE) {

 #' Get the number of genes for a \code{\linkS4class{GRN}} object.

 #' 

-#' Return the number of genes (all or only non-filtered ones) that are defined in the \code{\linkS4class{GRN}} object.

+#' Returns the number of genes (all or only non-filtered ones) from the provided RNA-seq data in the \code{\linkS4class{GRN}} object.

 #'

 #' @template GRN

 #' @param filter TRUE or FALSE. Default TRUE. Should genes marked as filtered be included in the count?

 #' @return Integer. Number of genes that are defined in the \code{\linkS4class{GRN}} object, either by excluding (filter = TRUE) or including (filter = FALSE) genes that are currently marked as \emph{filtered}.

+#' @seealso \code{\link{nTFs}}

+#' @seealso \code{\link{nPeaks}}

 #' @examples

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -288,7 +291,6 @@ nPeaks <- function(GRN, filter = TRUE) {

 #' nGenes(GRN, filter = FALSE)

 #' @export

 #' @aliases genes

-#' @rdname genes-methods

 nGenes <- function(GRN, filter = TRUE) {

   checkmate::assertClass(GRN, "GRN")

@@ -306,4 +308,32 @@ nGenes <- function(GRN, filter = TRUE) {

   nGenes

+}

+

+

+

+#' Get the number of TFs for a \code{\linkS4class{GRN}} object.

+#' 

+#' Returns the number of TFs  from the provided TFBS data in the \code{\linkS4class{GRN}} object.

+#'

+#' @template GRN

+#' @return Integer. Number of TFs that are defined in the \code{\linkS4class{GRN}} object.

+#' @seealso \code{\link{nGenes}}

+#' @seealso \code{\link{nPeaks}}

+#' @examples

+#' # See the Workflow vignette on the GRaNIE website for examples

+#' GRN = loadExampleObject()

+#' nTFs(GRN)

+#' @export

+#' @aliases TFs

+nTFs <- function(GRN) {

+    

+    checkmate::assertClass(GRN, "GRN")

+    

+    if (is.null(GRN@data$TFs$translationTable)) {

+        return(NA)

+    }

+    

+    nrow(GRN@data$TFs$translationTable)

+    

 }

\ No newline at end of file

@@ -1,24 +1,42 @@

 ######## Init GRN ########

-#' Initialize a \code{\linkS4class{GRN}} object

+#' Create and initialize a \code{\linkS4class{GRN}} object.

+#' 

+#' Executing this function is the very first step in the *GRaNIE* workflow. After its execution, data can be added to the object. 

+#' \strong{Depending on the genome assembly version, additional genome annotation packages are required, as follows:} 

+#' For \code{hg19} and \code{hg38},

+#' the packages \code{org.Hs.eg.db} as well as \code{BSgenome.Hsapiens.UCSC.hg19}+\code{TxDb.Hsapiens.UCSC.hg19.knownGene} or 

+#' \code{BSgenome.Hsapiens.UCSC.hg38}+\code{TxDb.Hsapiens.UCSC.hg38.knownGene} are required, respectively. 

+#' For \code{mm9} and \code{mm10},

+#' the packages \code{org.Mm.eg.db} as well as \code{BSgenome.Mmusculus.UCSC.mm9}+\code{TxDb.Mmusculus.UCSC.mm9.knownGene} or 

+#' \code{Mmusculus.UCSC.mm10}+\code{TxDb.Mmusculus.UCSC.mm10.knownGene} are required, respectively.

+#' For more information, see the error message if any of these packages is missing or the 

+#' \href{https://grp-zaugg.embl-community.io/GRaNIE/articles/GRaNIE_packageDetails.html#installation}{Package Details Vignette}.

 #' 

 #' @export

-#' @param objectMetadata List. Default \code{list()}. Optional (named) list with an arbitrary number of elements, all of which capture metadata for the object. This is mainly used to distinguish GRN objects from one another by storing object-specific metadata along with the data.

-#' @param outputFolder Output folder, either absolute or relative to the current working directory. No default. Default output folder where all pipeline output will be put unless specified otherwise. We recommend specifying an absolute path. Note that for Windows-based systems, the path must be correctly specified with "/" as path separator.

-#' @param genomeAssembly Character. No default. The genome assembly of all data that to be used within this object. Currently, supported genomes are: \code{hg19}, \code{hg38}, and \code{mm10}.

+#' @param objectMetadata List. Default \code{list()}. Optional (named) list with an arbitrary number of elements, all of which 

+#' capture metadata for the object. This is mainly used to distinguish GRN objects from one another by storing object-specific metadata along with the data.

+#' @param outputFolder Output folder, either absolute or relative to the current working directory. Default \code{"."}. 

+#' Default output folder where all pipeline output will be put unless specified otherwise. We recommend specifying an absolute path. 

+#' Note that for Windows-based systems, the path must be correctly specified with "/" as path separator.

+#' @param genomeAssembly Character. No default. The genome assembly of all data that to be used within this object. 

+#' Currently, supported genomes are: \code{hg19}, \code{hg38}, \code{mm9} and \code{mm10}. See function description for further information and notes.

 #' @return Empty \code{\linkS4class{GRN}} object

 #' @examples 

 #' meta.l = list(name = "exampleName", date = "01.03.22")

 #' GRN = initializeGRN(objectMetadata = meta.l, outputFolder = "output", genomeAssembly = "hg38")

 #' @export

-#' @import tidyverse

 #' @importFrom stats sd median cor cor.test quantile

 initializeGRN <- function(objectMetadata = list(),

-                          outputFolder, 

+                          outputFolder = ".", 

                           genomeAssembly) {

+  start = Sys.time()   

+    

   checkmate::assert(checkmate::checkNull(objectMetadata), checkmate::checkList(objectMetadata))

-  checkmate::assertSubset(genomeAssembly, c("hg19","hg38", "mm10"))

+  checkmate::assertChoice(genomeAssembly, c("hg19","hg38", "mm9", "mm10"))

+

+  .checkAndLoadPackagesGenomeAssembly(genomeAssembly)

   # Create the folder first if not yet existing

   if (!dir.exists(outputFolder)) {

@@ -51,8 +69,7 @@ initializeGRN <- function(objectMetadata = list(),

   packageName = utils::packageName()

   par.l$packageVersion = ifelse(is.null(packageName), NA, paste0(utils::packageVersion(packageName)[[1]], collapse = "."))

   par.l$genomeAssembly = genomeAssembly

-  

-  .checkAndLoadPackagesGenomeAssembly(genomeAssembly)

+

   # Make an internal subslot

   par.l$internal = list()

@@ -92,29 +109,47 @@ 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."))

+  .printExecutionTime(start, prefix = "")

   GRN

 }

 ######## Add and filter data ########

-#' Add data to a \code{\linkS4class{GRN}} object

+#' Add data to a \code{\linkS4class{GRN}} object.

+#' 

+#' This function adds both RNA and peak data to a \code{\linkS4class{GRN}} object, along with data normalization.

+#' In addition, and highly recommended, sample metadata can be optionally provided.

+#' 

+#' If the \code{ChIPseeker} package is installed, additional peak annotation is provided in the annotation slot and a peak annotation QC plot is produced as part of peak-gene QC.

+#' This is fully optional, however, and has no consequences for downstream functions.

 #' 

 #' @export

 #' @template GRN 

-#' @param counts_peaks Data frame. No default. Counts for the peaks, with raw or normalized counts for each peak (rows) across all samples (columns). In addition to the count data, it must also contain one ID column with a particular format, see the argument \code{idColumn_peaks} below. Row names are ignored, column names must be set to the sample names and must match those from the RNA counts and the sample metadata table.

-#' @param normalization_peaks Character. Default \code{DESeq_sizeFactor}. Normalization procedure for peak data. Must be one of \code{DESeq_sizeFactor}, \code{none}, or \code{quantile}.

-#' @param idColumn_peaks Character. Default \code{peakID}. Name of the column in the counts_peaks data frame that contains peak IDs. The required format must be {chr}:{start}-{end}", with {chr} denoting the abbreviated chromosome name, and {start} and {end} the begin and end of the peak coordinates, respectively. End must be bigger than start. Examples for valid peak IDs are \code{chr1:400-800} or \code{chrX:20-25}.

-#' @param counts_rna Data frame. No default. Counts for the RNA-seq data, with raw or normalized counts for each gene (rows) across all samples (columns). In addition to the count data, it must also contain one ID column with a particular format, see the argument \code{idColumn_rna} below. Row names are ignored, column names must be set to the sample names and must match those from the RNA counts and the sample metadata table.

-#' @param normalization_rna Character. Default \code{quantile}. Normalization procedure for peak data. Must be one of "DESeq_sizeFactor", "none", or "quantile"

+#' @param counts_peaks Data frame. No default. Counts for the peaks, with raw or normalized counts for each peak (rows) across all samples (columns). 

+#' In addition to the count data, it must also contain one ID column with a particular format, see the argument \code{idColumn_peaks} below. 

+#' Row names are ignored, column names must be set to the sample names and must match those from the RNA counts and the sample metadata table.

+#' @param normalization_peaks Character. Default \code{DESeq_sizeFactor}. 

+#' Normalization procedure for peak data. Must be one of \code{DESeq_sizeFactor}, \code{none}, or \code{quantile}.

+#' @param idColumn_peaks Character. Default \code{peakID}. Name of the column in the counts_peaks data frame that contains peak IDs. 

+#' The required format must be {chr}:{start}-{end}", with {chr} denoting the abbreviated chromosome name, and {start} and {end} the begin and end 

+#' of the peak coordinates, respectively. End must be bigger than start. Examples for valid peak IDs are \code{chr1:400-800} or \code{chrX:20-25}.

+#' @param counts_rna Data frame. No default. Counts for the RNA-seq data, with raw or normalized counts for each gene (rows) across all samples (columns). 

+#' In addition to the count data, it must also contain one ID column with a particular format, see the argument \code{idColumn_rna} below. 

+#' Row names are ignored, column names must be set to the sample names and must match those from the RNA counts and the sample metadata table.

+#' @param normalization_rna Character. Default \code{quantile}. Normalization procedure for peak data. 

+#' Must be one of "DESeq_sizeFactor", "none", or "quantile". For "quantile", \code{limma::normalizeQuantiles} is used for normalization.

 #' @param idColumn_RNA Character. Default \code{ENSEMBL}. Name of the column in the \code{counts_rna} data frame that contains Ensembl IDs.

-#' @param sampleMetadata Data frame. Default \code{NULL}. Optional, additional metadata for the samples, such as age, sex, gender etc. If provided, the @seealso [plotPCA_all()] function can then incorporate and plot it. Sample names must match with those from both peak and RNA-Seq data. The first column is expected to contain the sample IDs, the actual column name is irrelevant.

-#' @param allowOverlappingPeaks \code{TRUE} or \code{FALSE}. Default \code{FALSE}. Should overlapping peaks be allowed (then only a warning is issued when overlapping peaks are found) or (the default) should an error be raised?

+#' @param sampleMetadata Data frame. Default \code{NULL}. Optional, additional metadata for the samples, such as age, sex, gender etc. 

+#' If provided, the @seealso [plotPCA_all()] function can then incorporate and plot it. Sample names must match with those from both peak and RNA-Seq data. The first column is expected to contain the sample IDs, the actual column name is irrelevant.

+#' @param allowOverlappingPeaks \code{TRUE} or \code{FALSE}. Default \code{FALSE}. Should overlapping peaks be allowed (then only a warning is issued 

+#' when overlapping peaks are found) or (the default) should an error be raised?

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.

+#' @return An updated \code{\linkS4class{GRN}} object, with added data from this function (e.g., slots \code{GRN@data$peaks} and \code{GRN@data$RNA})

+#' @seealso \code{\link{plotPCA_all}}

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

-#' # library(tidyverse)

+#' # library(readr)

 #' # rna.df   = read_tsv("https://www.embl.de/download/zaugg/GRaNIE/rna.tsv.gz")

 #' # peaks.df = read_tsv("https://www.embl.de/download/zaugg/GRaNIE/peaks.tsv.gz")

 #' # meta.df  = read_tsv("https://www.embl.de/download/zaugg/GRaNIE/sampleMetadata.tsv.gz")

@@ -127,15 +162,18 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

                     allowOverlappingPeaks= FALSE,

                     forceRerun = FALSE) {

-  GRN = .addFunctionLogToObject(GRN)      

+  start = Sys.time()

   checkmate::assertClass(GRN, "GRN")

+  GRN = .addFunctionLogToObject(GRN)      

+  

   checkmate::assertDataFrame(counts_peaks, min.rows = 1, min.cols = 2)

   checkmate::assertDataFrame(counts_rna, min.rows = 1, min.cols = 2)

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

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

   checkmate::assertChoice(normalization_peaks, c("none", "DESeq_sizeFactor", "quantile"))

   checkmate::assertChoice(normalization_rna, c("none", "DESeq_sizeFactor", "quantile"))  

+  checkmate::assertFlag(allowOverlappingPeaks)

   checkmate::assertFlag(forceRerun)

   if (is.null(GRN@data$peaks$counts_norm) |

@@ -328,6 +366,8 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

   # Add gene annotation once

   GRN = .populateGeneAnnotation(GRN)

+  .printExecutionTime(start, prefix = "")

+  

   GRN

 }

@@ -369,10 +409,10 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

   paste0(ids_chr, ":", format(as.integer(start), scientific = FALSE), "-", format(as.integer(end), scientific= FALSE))

 }

-

+#' @importFrom biomaRt useEnsembl getBM

 .retrieveAnnotationData <- function(genomeAssembly) {

-    futile.logger::flog.info(paste0("Retrieving genome annotation data from biomaRt... This may take a while"))

+    futile.logger::flog.info(paste0("Retrieving genome annotation data from biomaRt for ", genomeAssembly, "... This may take a while"))

     params.l = .getBiomartParameters(genomeAssembly)

@@ -407,7 +447,9 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

         dplyr::mutate(chromosome_name = paste0("chr", chromosome_name)) %>%

         dplyr::rename(gene.chr = chromosome_name, gene.start = start_position, gene.end = end_position, 

                       gene.strand = strand, gene.ENSEMBL = ensembl_gene_id, gene.type = gene_biotype, gene.name = external_gene_name) %>%

+        tidyr::replace_na(list(gene.type = "unknown")) %>%

         dplyr::mutate_if(is.character, as.factor) %>%

+        dplyr::mutate(gene.type = dplyr::recode_factor(gene.type, lncRNA = "lincRNA")) %>%  # there seems to be a name change from lincRNA -> lncRNA, lets change it here 

         dplyr::mutate(gene.strand = factor(gene.strand, levels = c(1,-1), labels = c("+", "-")))

 }

@@ -458,11 +500,11 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

                                   peak.CV = CV_peaks)

   GRN@annotation$consensusPeaks = metadata_peaks

-  

+

   if (!is.installed("ChIPseeker")) {

-    message = paste0("The package ChIPseeker is currently not installed, which is needed for additional peak annotation that can be useful for further downstream analyses. ", 

-                     " You may want to install it and re-run this function. However, this is optional and except for some missing additional annotation columns, there is no limitation.")

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

+      packageMessage = paste0("The package ChIPseeker is currently not installed, which is needed for additional peak annotation that can be useful for further downstream analyses. ", 

+                              " You may want to install it and re-run this function. However, this is optional and except for some missing additional annotation columns, there is no limitation.")

+      .checkPackageInstallation("ChIPseeker", packageMessage, isWarning = TRUE)

   } else {

     futile.logger::flog.info(paste0(" Retrieve peak annotation using ChipSeeker. This may take a while"))

@@ -545,7 +587,8 @@ 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::full_join(genomeAnnotation.df, by = c("gene.ENSEMBL")) %>%

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

   GRN@annotation$genes = metadata_rna

@@ -560,7 +603,6 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 }

-#' @importFrom IRanges width

 #' @importFrom rlang .data `:=`

 .calcGCContentPeaks <- function(GRN) {

@@ -582,7 +624,7 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

   GC_content.df = GC_content.df %>%

     tibble::as_tibble() %>%

-    dplyr::mutate(length = IRanges::width(query),

+    dplyr::mutate(length = GenomicRanges::width(query),

                   peak.ID = query$peakID,

                   GC_class = cut(`G|C`, breaks = seq(0,1,0.1), include.lowest = TRUE, ordered_result = TRUE))

@@ -611,22 +653,22 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

   GRN

 }

-#' Filter data from a \code{\linkS4class{GRN}} object

+#' Filter RNA-seq and/or peak data from a \code{\linkS4class{GRN}} object

 #' 

 #' @template GRN 

-#' @param minNormalizedMean_peaks Numeric or \code{NULL}. Default 5. Minimum mean across all samples for a peak to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

-#' @param maxNormalizedMean_peaks Numeric or \code{NULL}. Default \code{NULL}. Maximum mean across all samples for a peak to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

-#' @param minNormalizedMeanRNA Numeric or \code{NULL}. Default 5. Minimum mean across all samples for a gene to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

-#' @param maxNormalizedMeanRNA Numeric or \code{NULL}. Default \code{NULL}. Maximum mean across all samples for a gene to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

+#' @param minNormalizedMean_peaks Numeric[0,] or \code{NULL}. Default 5. Minimum mean across all samples for a peak to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

+#' @param maxNormalizedMean_peaks Numeric[0,] or \code{NULL}. Default \code{NULL}. Maximum mean across all samples for a peak to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

+#' @param minNormalizedMeanRNA Numeric[0,] or \code{NULL}. Default 5. Minimum mean across all samples for a gene to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

+#' @param maxNormalizedMeanRNA Numeric[0,] or \code{NULL}. Default \code{NULL}. Maximum mean across all samples for a gene to be retained for the normalized counts table. Set to \code{NULL} for not applying the filter.

 #' @param chrToKeep_peaks Character vector or \code{NULL}. Default \code{NULL}. Vector of chromosomes that peaks are allowed to come from. This filter can be used to filter sex chromosomes from the peaks, for example (e.g, \code{c(paste0("chr", 1:22), "chrX", "chrY")})

-#' @param minSize_peaks Integer or \code{NULL}. Default \code{NULL}. Minimum peak size (width, end - start) for a peak to be retained. Set to \code{NULL} for not applying the filter.

-#' @param maxSize_peaks Integer or \code{NULL}. Default 10000. Maximum peak size (width, end - start) for a peak to be retained. Set to \code{NULL} for not applying the filter.

-#' @param minCV_peaks Numeric or \code{NULL}. Default \code{NULL}. Minimum CV (coefficient of variation, a unitless measure of variation) for a peak to be retained. Set to \code{NULL} for not applying the filter.

-#' @param maxCV_peaks Numeric or \code{NULL}. Default \code{NULL}. Maximum CV (coefficient of variation, a unitless measure of variation) for a peak to be retained. Set to \code{NULL} for not applying the filter.

-#' @param minCV_genes Numeric or \code{NULL}. Default \code{NULL}. Minimum CV (coefficient of variation, a unitless measure of variation) for a gene to be retained. Set to \code{NULL} for not applying the filter.

-#' @param maxCV_genes Numeric or \code{NULL}. Default \code{NULL}. Maximum CV (coefficient of variation, a unitless measure of variation) for a gene to be retained. Set to \code{NULL} for not applying the filter.

+#' @param minSize_peaks Integer[1,] or \code{NULL}. Default \code{NULL}. Minimum peak size (width, end - start) for a peak to be retained. Set to \code{NULL} for not applying the filter.

+#' @param maxSize_peaks Integer[1,] or \code{NULL}. Default 10000. Maximum peak size (width, end - start) for a peak to be retained. Set to \code{NULL} for not applying the filter.

+#' @param minCV_peaks Numeric[0,] or \code{NULL}. Default \code{NULL}. Minimum CV (coefficient of variation, a unitless measure of variation) for a peak to be retained. Set to \code{NULL} for not applying the filter.

+#' @param maxCV_peaks Numeric[0,] or \code{NULL}. Default \code{NULL}. Maximum CV (coefficient of variation, a unitless measure of variation) for a peak to be retained. Set to \code{NULL} for not applying the filter.

+#' @param minCV_genes Numeric[0,] or \code{NULL}. Default \code{NULL}. Minimum CV (coefficient of variation, a unitless measure of variation) for a gene to be retained. Set to \code{NULL} for not applying the filter.

+#' @param maxCV_genes Numeric[0,] or \code{NULL}. Default \code{NULL}. Maximum CV (coefficient of variation, a unitless measure of variation) for a gene to be retained. Set to \code{NULL} for not applying the filter.

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.

+#' @return An updated \code{\linkS4class{GRN}} object, with added data from this function.

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -640,9 +682,12 @@ filterData <- function (GRN,

                         minCV_peaks = NULL, maxCV_peaks = NULL,

                         minCV_genes = NULL, maxCV_genes = NULL,

                         forceRerun = FALSE) {

-  GRN = .addFunctionLogToObject(GRN) 

+  start = Sys.time()

+    

   checkmate::assertClass(GRN, "GRN")

+  GRN = .addFunctionLogToObject(GRN) 

+

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

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

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

@@ -717,6 +762,8 @@ filterData <- function (GRN,

     GRN@config$isFiltered = TRUE

   } 

+  .printExecutionTime(start, prefix = "")

+  

   GRN

 }

@@ -887,17 +934,18 @@ filterData <- function (GRN,

 #' @param TFs Character vector. Default \code{all}. Vector of TF names to include. The special keyword \code{all} can be used to include all TF found in the folder as specified by \code{motifFolder}. If \code{all} is specified anywhere, all TFs will be included. TF names must otherwise match the file names that are found in the folder, without the file suffix.

 #' @param filesTFBSPattern Character. Default \code{"_TFBS"}. Suffix for the file names in the TFBS folder that is not part of the TF name. Can be empty. For example, for the TF CTCF, if the file is called \code{CTCF.all.TFBS.bed}, set this parameter to \code{".all.TFBS"}.

 #' @param fileEnding Character. Default \code{".bed"}. File ending for the files from the motif folder.

-#' @param nTFMax \code{NULL} or integer. Default \code{NULL}. Maximal number of TFs to import. Can be used for testing purposes, e.g., setting to 5 only imports 5 TFs even though the whole \code{motifFolder} has many more TFs defined.

+#' @param nTFMax \code{NULL} or Integer[1,]. Default \code{NULL}. Maximal number of TFs to import. Can be used for testing purposes, e.g., setting to 5 only imports 5 TFs even though the whole \code{motifFolder} has many more TFs defined.

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.

+#' @return An updated \code{\linkS4class{GRN}} object, with additional information added from this function (\code{GRN@data$TFs$translationTable} in particular)

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' @export

 addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPattern = "_TFBS", fileEnding = ".bed", forceRerun = FALSE) {

-  GRN = .addFunctionLogToObject(GRN)

-  

+  start = Sys.time()

   checkmate::assertClass(GRN, "GRN")

+  GRN = .addFunctionLogToObject(GRN)

+

   checkmate::assertDirectoryExists(motifFolder)

   checkmate::assertCharacter(TFs, min.len = 1)

   checkmate::assert(checkmate::testNull(nTFMax), checkmate::testIntegerish(nTFMax, lower = 1))

@@ -927,6 +975,8 @@ addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPatte

   } 

+  .printExecutionTime(start, prefix = "")

+  

   GRN

 }

@@ -1006,7 +1056,7 @@ addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPatte

 #' @template GRN

 #' @template nCores

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function. 

+#' @return An updated \code{\linkS4class{GRN}} object, with added data from this function (\code{GRN@data$TFs$TF_peak_overlap} in particular)

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -1014,9 +1064,11 @@ addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPatte

 #' @export

 overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

+  start = Sys.time()

+    

+  checkmate::assertClass(GRN, "GRN")

   GRN = .addFunctionLogToObject(GRN)

-  checkmate::assertClass(GRN, "GRN")

   checkmate::assertIntegerish(nCores, lower = 1)

   checkmate::assertFlag(forceRerun)

@@ -1092,6 +1144,8 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

   } 

+  .printExecutionTime(start, prefix = "")

+  

   GRN

 }

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

     dplyr::group_by(peakID) %>%

     dplyr::slice(which.min(distance)) %>%

     #arrange(distance, .by_group = TRUE) %>%

-    # top_n(n = 2, desc(distance)) %>%

+    # top_n(n = 2, dplyr::desc(distance)) %>%

     dplyr::ungroup()

   futile.logger::flog.info(paste0(" Calculating intersection for TF ", TFCur, " finished. Number of overlapping TFBS after filtering: ", nrow(final.df)))

@@ -1226,7 +1280,7 @@ 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(ENSEMBL %in% colnames(sort.cor.m))

+  HOCOMOCO_mapping.exp.filt = HOCOMOCO_mapping.exp %>% dplyr::filter(.data$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)

@@ -1239,8 +1293,8 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

   peak_TF_overlapCur.df = .asMatrixFromSparse(GRN@data$TFs$TF_peak_overlap, convertZero_to_NA = FALSE) %>% 

     tibble::as_tibble() %>%

-    dplyr::filter(!isFiltered) %>% 

-    dplyr::select(-isFiltered) 

+    dplyr::filter(!.data$isFiltered) %>% 

+    dplyr::select(-.data$isFiltered) 

   if (perm > 0 & shuffle) {

     peak_TF_overlapCur.df = .shuffleRowsPerColumn(peak_TF_overlapCur.df)

@@ -1257,18 +1311,22 @@ overlapPeaksAndTFBS <- function(GRN, nCores = 2, forceRerun = FALSE) {

 #' Add TF activity data to GRN object using a simplified procedure for estimating it. EXPERIMENTAL.

 #' 

+#' We do not yet provide full support for this function. It is currently being tested. Use at our own risk.

+#' 

 #' @template GRN

 #' @template normalization_TFActivity

 #' @template name_TFActivity

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function. 

+#' @return An updated \code{\linkS4class{GRN}} object, with added data from this function

+#' (\code{GRN@data$TFs[[name]]} in particular, with \code{name} referring to the value of tje \code{name} parameter) 

 addData_TFActivity <- function(GRN, normalization = "cyclicLoess", name = "TF_activity", forceRerun = FALSE) {

+  checkmate::assertClass(GRN, "GRN")

   GRN = .addFunctionLogToObject(GRN)

   start = Sys.time()

-  checkmate::assertClass(GRN, "GRN")

+  

   checkmate::assertChoice(normalization, c("cyclicLoess", "sizeFactors", "quantile", "none"))

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

   checkmate::assertFlag(forceRerun)

@@ -1365,6 +1423,11 @@ addData_TFActivity <- function(GRN, normalization = "cyclicLoess", name = "TF_ac

   if (normalization == "cyclicLoess") {

     futile.logger::flog.info(paste0(" Normalizing data using cyclic LOESS"))

+      

+      

+    packageMessage = paste0("The package csaw is not installed but required for the cyclic LOESS normalization. Please install it and re-run this function or change the normalization type (if possible).")

+    .checkPackageInstallation("csaw", packageMessage)   

+      

     # Perform a cyclic loess normalization

     # We use a slighlty more complicated setup to derive size factors for library normalization

     # Instead of just determining the size factors in DeSeq2 via cirtual samples, we use 

@@ -1413,6 +1476,8 @@ addData_TFActivity <- function(GRN, normalization = "cyclicLoess", name = "TF_ac

 #' Import externally derived TF Activity data. EXPERIMENTAL.

 #' 

+#' We do not yet provide full support for this function. It is currently being tested. Use at our own risk.

+#' 

 #' @template GRN

 #' @param data Data frame. No default. Data with TF data.

 #' @template name_TFActivity

@@ -1420,18 +1485,19 @@ addData_TFActivity <- function(GRN, normalization = "cyclicLoess", name = "TF_ac

 #' @param nameColumn Character. Default \code{TF.name}. Must be unique for each TF / row.

 #' @template normalization_TFActivity

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.  

+#' @return An updated \code{\linkS4class{GRN}} object, with added data from this function.  

 importTFData <- function(GRN, data, name, idColumn = "ENSEMBL", nameColumn = "TF.name", normalization = "none", forceRerun = FALSE) {

+  checkmate::assertClass(GRN, "GRN")

   GRN = .addFunctionLogToObject(GRN)

   start = Sys.time()

-  checkmate::assertClass(GRN, "GRN")

+  

   checkmate::assertDataFrame(data, min.cols = 2, min.rows = 1)

-  checkmate::assertSubset(idColumn, colnames(data))

-  checkmate::assertSubset(nameColumn, colnames(data))

-  checkmate::assertSubset(GRN@config$sharedSamples, colnames(data))

+  checkmate::assertChoice(idColumn, colnames(data))

+  checkmate::assertChoice(nameColumn, colnames(data))

+  checkmate::assertSubset(GRN@config$sharedSamples, colnames(data), empty.ok = FALSE)

   checkmate::assertCharacter(name, min.chars = 1, any.missing = FALSE, len = 1)

   checkmate::assertChoice(normalization, c("cyclicLoess", "sizeFactors", "quantile", "none"))

@@ -1482,7 +1548,7 @@ importTFData <- function(GRN, data, name, idColumn = "ENSEMBL", nameColumn = "TF

     countsNorm$ENSEMBL = data$ENSEMBL

     nRowBefore = nrow(countsNorm)

-    countsNorm.subset = dplyr::filter(countsNorm, ENSEMBL %in% GRN@data$TFs$translationTable$ENSEMBL)

+    countsNorm.subset = dplyr::filter(countsNorm, .data$ENSEMBL %in% GRN@data$TFs$translationTable$ENSEMBL)

     nRowAfter = nrow(countsNorm.subset)

     if (nRowAfter == 0) {

       message = "No rows overlapping with translation table, check ENSEMBL IDs."

@@ -1535,14 +1601,14 @@ importTFData <- function(GRN, data, name, idColumn = "ENSEMBL", nameColumn = "TF

 #' Run the activator-repressor classification for the TFs for a \code{\linkS4class{GRN}} object

 #' 

 #' @template GRN

-#' @param significanceThreshold_Wilcoxon Numeric between 0 and 1. Default 0.05. Significance threshold for Wilcoxon test that is run in the end for the final classification. See the Vignette and *diffTF* paper for details.

-#' @param plot_minNoTFBS_heatmap Integer. Default 100. Minimum number of TFBS for a TF to be included in the heatmap that is part of the output of this function.

+#' @param significanceThreshold_Wilcoxon Numeric[0,1]. Default 0.05. Significance threshold for Wilcoxon test that is run in the end for the final classification. See the Vignette and *diffTF* paper for details.

+#' @param plot_minNoTFBS_heatmap Integer[1,]. Default 100. Minimum number of TFBS for a TF to be included in the heatmap that is part of the output of this function.

 #' @param deleteIntermediateData \code{TRUE} or \code{FALSE}.  Default \code{TRUE}. Should intermediate data be deleted before returning the object after a successful run? Due to the size of the produced intermediate data, we recommend setting this to \code{TRUE}, but if memory or object size are not an issue, the information can also be kept.

 #' @template plotDiagnosticPlots

 #' @template outputFolder

 #' @template corMethod

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.  

+#' @return An updated \code{\linkS4class{GRN}} object, with additional information added from this function. 

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' # GRN = loadExampleObject()

@@ -1554,9 +1620,12 @@ AR_classification_wrapper<- function (GRN, significanceThreshold_Wilcoxon = 0.05

                                       corMethod = "pearson",

                                       forceRerun = FALSE) {

+  start = Sys.time()

+    

+  checkmate::assertClass(GRN, "GRN")

   GRN = .addFunctionLogToObject(GRN)

-  checkmate::assertClass(GRN, "GRN")

+  

   checkmate::assertNumber(significanceThreshold_Wilcoxon, lower = 0, upper = 1)

   checkmate::assertNumber(plot_minNoTFBS_heatmap, lower = 1)

   checkmate::assertFlag(deleteIntermediateData)

@@ -1730,6 +1799,7 @@ AR_classification_wrapper<- function (GRN, significanceThreshold_Wilcoxon = 0.05

   } # end of for each connection type

+  .printExecutionTime(start, prefix = "")

   GRN

@@ -1753,6 +1823,8 @@ AR_classification_wrapper<- function (GRN, significanceThreshold_Wilcoxon = 0.05

 #' Add TF-peak connections to a \code{\linkS4class{GRN}} object

 #' 

+#' After the execution of this function, QC plots can be plotted with the function \code{\link{plotDiagnosticPlots_TFPeaks}} unless this has already been done by default due to \code{plotDiagnosticPlots = TRUE}

+#' 

 #' @template GRN 

 #' @template plotDiagnosticPlots

 #' @template plotDetails

@@ -1760,13 +1832,14 @@ AR_classification_wrapper<- function (GRN, significanceThreshold_Wilcoxon = 0.05

 #' @template corMethod

 #' @param connectionTypes Character vector. Default \code{expression}. Vector of connection types to include for the TF-peak connections. If an additional connection type is specified here, it has to be available already within the object (EXPERIMENTAL). See the function \code{\link{addData_TFActivity}} for details.

 #' @param removeNegativeCorrelation  Vector of \code{TRUE} or \code{FALSE}. Default \code{FALSE}. EXPERIMENTAL. Must be a logical vector of the same length as the parameter \code{connectionType}. Should negatively correlated TF-peak connections be removed for the specific connection type? For connection type expression, the default is \code{FALSE}, while for any TF Activity related connection type, we recommend setting this to \code{TRUE}.  

-#' @param maxFDRToStore Numeric. Default 0.3. Maximum TF-peak FDR value to permanently store a particular TF-peak connection in the object? This parameter has a large influence on the overall memory size of the object, and we recommend not storing connections with a high FDR due to their sheer number.

+#' @param maxFDRToStore Numeric[0,1]. Default 0.3. Maximum TF-peak FDR value to permanently store a particular TF-peak connection in the object? This parameter has a large influence on the overall memory size of the object, and we recommend not storing connections with a high FDR due to their sheer number.

 #' @param addForPermuted \code{TRUE} or \code{FALSE}.  Default \code{FALSE}. Add connections also for permuted data. Leave at \code{TRUE} unless you know what you are doing.

 #' @param useGCCorrection \code{TRUE} or \code{FALSE}.  Default \code{FALSE}. EXPERIMENTAL. Should a GC-matched background be used when calculating FDRs?

-#' @param percBackground_size Numeric (0 to 100). Default 75. EXPERIMENTAL. Description will follow. Only relevant if \code{useGCCorrection} is set to \code{TRUE}, ignored otherwise.

+#' @param percBackground_size Numeric[0,100]. Default 75. EXPERIMENTAL. Description will follow. Only relevant if \code{useGCCorrection} is set to \code{TRUE}, ignored otherwise.

 #' @param percBackground_resample \code{TRUE} or \code{FALSE}.  Default \code{TRUE}. EXPERIMENTAL. Should resampling be enabled for those GC bins for which not enough background peaks are available?. Only relevant if \code{useGCCorrection} is set to \code{TRUE}, ignored otherwise.

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function.

+#' @seealso \code{\link{plotDiagnosticPlots_TFPeaks}}

+#' @return An updated \code{\linkS4class{GRN}} object, with additional information added from this function. 

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -1781,16 +1854,18 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

                                     useGCCorrection = FALSE, percBackground_size = 75, percBackground_resample = TRUE,

                                     forceRerun = FALSE) {

-  GRN = .addFunctionLogToObject(GRN)

-  

+  start = Sys.time()

+

   checkmate::assertClass(GRN, "GRN")

+  GRN = .addFunctionLogToObject(GRN)

+

   checkmate::assertFlag(plotDiagnosticPlots)

   checkmate::assertFlag(plotDetails)

   checkmate::assertChoice(corMethod, c("pearson", "spearman"))

   checkmate::assertFlag(addForPermuted)

   GRN = .checkAndUpdateConnectionTypes(GRN) # For compatibility with older versions

-  checkmate::assertSubset(connectionTypes, GRN@config$TF_peak_connectionTypes)

+  checkmate::assertSubset(connectionTypes, GRN@config$TF_peak_connectionTypes, empty.ok = FALSE)

   checkmate::assertLogical(removeNegativeCorrelation, any.missing = FALSE, len = length(connectionTypes))

@@ -1852,6 +1927,8 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

   }

+  .printExecutionTime(start, prefix = "")

+  

   GRN

 }

@@ -1928,7 +2005,7 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

     pb <- progress::progress_bar$new(total = length(allTF))

-    peaksFiltered = GRN@data$peaks$consensusPeaks %>% dplyr::filter(!isFiltered) 

+    peaksFiltered = GRN@data$peaks$consensusPeaks %>% dplyr::filter(!.data$isFiltered) 

     if (!useGCCorrection) {

       minPerc = 100

@@ -2036,7 +2113,7 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

         peakIDsSel = c()

         for (i in seq_len(nrow(GC_classes_foreground.df))) {

-          peaksBackgroundGCCur =  peaksBackground %>% dplyr::filter(GC_class == GC_classes_foreground.df$GC_class[i])

+          peaksBackgroundGCCur =  peaksBackground %>% dplyr::filter(.data$GC_class == GC_classes_foreground.df$GC_class[i])

           if (nrow( peaksBackgroundGCCur) == 0) {

             next

@@ -2180,7 +2257,7 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

         if (connectionTypeCur %in% connectionTypes_removeNegCor ) {

           # futile.logger::flog.info(paste0(" Remove negatively correlated TF-peak pairs for connection type ", connectionTypeCur))

-          cor.peak.tf = dplyr::filter(cor.peak.tf, TF_peak.r >= 0)

+          cor.peak.tf = dplyr::filter(cor.peak.tf, .data$TF_peak.r >= 0)

         }

@@ -2218,10 +2295,10 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

         # DISCARD other rows altogether

         # Left join here is what we want, as we need this df only for "real" data

         tblFilt.df = dplyr::left_join(cor.peak.tf, fdr.curve, by = "TF_peak.r_bin") %>%

-          dplyr::filter(TF_peak.fdr <= maxFDRToStore) %>%

-          dplyr::select(TF.name, TF_peak.r_bin, TF_peak.r, TF_peak.fdr, 

-                        TF_peak.fdr_orig, peak.ID, TF_peak.fdr_direction, 

-                        TF_peak.connectionType, tidyselect::contains("value"))

+          dplyr::filter(.data$TF_peak.fdr <= maxFDRToStore) %>%

+          dplyr::select(.data$TF.name, .data$TF_peak.r_bin, .data$TF_peak.r, .data$TF_peak.fdr, 

+                        .data$TF_peak.fdr_orig, .data$peak.ID, .data$TF_peak.fdr_direction, 

+                        .data$TF_peak.connectionType, tidyselect::contains("value"))

         if (!plotDetails) {

@@ -2277,7 +2354,7 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

       requiredNoPeaks = round(n_rel * targetNoPeaks, 0)

       # Check in background

       availableNoPeaks = GC_classes_background.df %>% 

-        dplyr::filter(GC_class == GC_class.cur) %>%

+        dplyr::filter(.data$GC_class == GC_class.cur) %>%

         dplyr::pull(.data$n)

       #futile.logger::flog.info(paste0(" GC.class ", GC.class.cur, ": Required: ", requiredNoPeaks, ", available: ", availableNoPeaks))

       if ( availableNoPeaks < requiredNoPeaks) {

@@ -2309,6 +2386,8 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

 #' Add peak-gene connections to a \code{\linkS4class{GRN}} object

 #' 

+#' After the execution of this function, QC plots can be plotted with the function \code{\link{plotDiagnosticPlots_peakGene}} unless this has already been done by default due to \code{plotDiagnosticPlots = TRUE}

+#' 

 #' @export

 #' @template GRN

 #' @param  overlapTypeGene Character. \code{"TSS"} or \code{"full"}. Default \code{"TSS"}. If set to \code{"TSS"}, only the TSS of the gene is used as reference for finding genes in the neighborhood of a peak. If set to \code{"full"}, the whole annotated gene (including all exons and introns) is used instead. 

@@ -2317,11 +2396,12 @@ addConnections_TF_peak <- function (GRN, plotDiagnosticPlots = TRUE, plotDetails

 #' @param TADs Data frame with TAD domains. Default \code{NULL}. If provided, the neighborhood of a peak is defined by the TAD domain the peak is in rather than a fixed-sized neighborhood. The expected format is a BED-like data frame with at least 3 columns in this particular order: chromosome, start, end, the 4th column is optional and will be taken as ID column. All additional columns as well as column names are ignored. For the first 3 columns, the type is checked as part of a data integrity check.

 #' @template nCores

 #' @template plotDiagnosticPlots

-#' @param plotGeneTypes List of character vectors. Default \code{list(c("all"), c("protein_coding"), c("protein_coding", "lincRNA"))}. Each list element may consist of one or multiple gene types that are plotted collectively in one PDF. The special keyword \code{"all"} denotes all gene types that are found (be aware: this typically contains 20+ gene types, see \url{https://www.gencodegenes.org/pages/biotypes.html} for details).

+#' @param plotGeneTypes List of character vectors. Default \code{list(c("all"), c("protein_coding"))}. Each list element may consist of one or multiple gene types that are plotted collectively in one PDF. The special keyword \code{"all"} denotes all gene types that are found (be aware: this typically contains 20+ gene types, see \url{https://www.gencodegenes.org/pages/biotypes.html} for details).

 #' @template outputFolder

 #' @template addRobustRegression

 #' @template forceRerun

-#' @return The same \code{\linkS4class{GRN}} object, with added data from this function in different flavors.

+#' @seealso \code{\link{plotDiagnosticPlots_peakGene}}

+#' @return An updated \code{\linkS4class{GRN}} object, with additional information added from this function. 

 #' @examples 

 #' # See the Workflow vignette on the GRaNIE website for examples

 #' GRN = loadExampleObject()

@@ -2330,25 +2410,36 @@ addConnections_peak_gene <- function(GRN, overlapTypeGene = "TSS", corMethod = "

                                      promoterRange = 250000, TADs = NULL,

                                      nCores = 4, 

                                      plotDiagnosticPlots = TRUE, 

-                                     plotGeneTypes = list(c("all"), c("protein_coding"), c("protein_coding", "lincRNA")), 

+                                     plotGeneTypes = list(c("all"), c("protein_coding")), 

                                      outputFolder = NULL,

                                      addRobustRegression = FALSE,

                                      forceRerun = FALSE) {

+  start = Sys.time() 

+    

+  checkmate::assertClass(GRN, "GRN")

   GRN = .addFunctionLogToObject(GRN)

-  checkmate::assertClass(GRN, "GRN")

+  

   checkmate::assertChoice(overlapTypeGene, c("TSS", "full"))

   checkmate::assertChoice(corMethod, c("pearson", "spearman"))

   checkmate::assertIntegerish(promoterRange, lower = 0)

   checkmate::assert(checkmate::testNull(TADs), checkmate::testDataFrame(TADs))