@@ -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: 0.99.0

+Version: 0.99.1

 Encoding: UTF-8

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

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

@@ -11,6 +11,7 @@ export(calculateCommunitiesEnrichment)

 export(calculateCommunitiesStats)

 export(calculateGeneralEnrichment)

 export(calculateTFEnrichment)

+export(changeOutputDirectory)

 export(deleteIntermediateData)

 export(filterData)

 export(filterGRNAndConnectGenes)

@@ -16,10 +16,10 @@

 initializeGRN <- function(objectMetadata = list(),

                           outputFolder, 

                           genomeAssembly) {

-

+  

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

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

-

+  

   # Create the folder first if not yet existing

   if (!dir.exists(outputFolder)) {

     dir.create(outputFolder)

@@ -44,7 +44,7 @@ initializeGRN <- function(objectMetadata = list(),

   GRN = .addFunctionLogToObject(GRN)

   GRN@config$isFiltered = FALSE

-

+  

   par.l = list()

   packageName = utils::packageName()

@@ -80,7 +80,7 @@ initializeGRN <- function(objectMetadata = list(),

   GRN@config$directories$outputRoot         = outputFolder

   GRN@config$directories$output_plots       = dir_output_plots 

   GRN@config$files$output_log               = paste0(outputFolder, "GRN.log")

-

+  

   .testExistanceAndCreateDirectoriesRecursively(c(outputFolder, dir_output_plots))

   checkmate::assertDirectory(outputFolder, access = "w")

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

                     counts_rna, normalization_rna = "quantile", idColumn_RNA = "ENSEMBL", sampleMetadata = NULL,

                     allowOverlappingPeaks= FALSE,

                     forceRerun = FALSE) {

-

+  

   GRN = .addFunctionLogToObject(GRN)      

   checkmate::assertClass(GRN, "GRN")

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

       forceRerun) {

     # Store raw peaks counts as DESeq object

-

+    

     #GRN@data$peaks$counts_raw = counts_peaks %>% dplyr::mutate(isFiltered = FALSE)

     # Normalize ID column names

@@ -152,31 +152,31 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     if (idColumn_RNA != "ENSEMBL") {

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

     }

-      

+    

     # Check ID columns for missing values and remove

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

     if (length(rna_missing_ID) > 0) {

-        message = paste0(" Found ", length(rna_missing_ID), " missing IDs in the ID column of the RNA counts. These rows will be removed.")

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

-        counts_rna = dplyr::slice(counts_rna, -rna_missing_ID)

+      message = paste0(" Found ", length(rna_missing_ID), " missing IDs in the ID column of the RNA counts. These rows will be removed.")

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

+      counts_rna = dplyr::slice(counts_rna, -rna_missing_ID)

     }

     peaks_missing_ID =  which(is.na(counts_peaks$peakID))

     if (length(peaks_missing_ID) > 0) {

-        message = paste0(" Found ", length(peaks_missing_ID), " missing IDs in the ID column of the peaks counts. These rows will be removed.")

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

-        counts_peaks = dplyr::slice(counts_peaks, -peaks_missing_ID)

+      message = paste0(" Found ", length(peaks_missing_ID), " missing IDs in the ID column of the peaks counts. These rows will be removed.")

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

+      counts_peaks = dplyr::slice(counts_peaks, -peaks_missing_ID)

     }

-      

+    

     # Remove potential scientific notation from peak IDs

     peaks_eNotation = which(grepl("e+", counts_peaks$peakID))

     if (length(peaks_eNotation) > 0) {

-        message = paste0("Found at least one peak (", paste0(counts_peaks$peakID[peaks_eNotation], collapse = ",") , ") for which the position contains the scientific notation, attempting to fix.")

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

-        counts_peaks$peakID[peaks_eNotation] = .removeScientificNotation_positions(counts_peaks$peakID[peaks_eNotation])

-   

-        

+      message = paste0("Found at least one peak (", paste0(counts_peaks$peakID[peaks_eNotation], collapse = ",") , ") for which the position contains the scientific notation, attempting to fix.")

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

+      counts_peaks$peakID[peaks_eNotation] = .removeScientificNotation_positions(counts_peaks$peakID[peaks_eNotation])

+      

+      

     }

     # Clean Ensembl IDs

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

       counts_rna = counts_rna %>%

         dplyr::group_by(ENSEMBL) %>%

         dplyr::summarise_if(is.numeric, sum) 

-        # dplyr::summarise_if(is.numeric, sum, .groups = 'drop') # the .drop caused an error with dplyr 1.0.5

+      # dplyr::summarise_if(is.numeric, sum, .groups = 'drop') # the .drop caused an error with dplyr 1.0.5

     }

-

+    

     # Normalize counts

     countsPeaks.norm.df  = .normalizeCounts(counts_peaks, method = normalization_peaks, idColumn = "peakID")

     countsRNA.norm.df  = .normalizeCounts(counts_rna, method = normalization_rna, idColumn = "ENSEMBL")

@@ -219,20 +219,20 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     # Create permutations for RNA

     futile.logger::flog.info(paste0( "Produce ", .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)) {

-        

+      

       futile.logger::flog.info("Parsing provided metadata...")

       GRN@data$metadata = sampleMetadata %>% dplyr::distinct() %>% tibble::set_tidy_names(syntactic = TRUE, quiet = TRUE)

       # Force the first column to be the ID column

       if ("sampleID" %in% colnames(GRN@data$metadata)) {

-          futile.logger::flog.warn("Renaming the first column to sampleID. However, this column already exists, it will be renamed accordingly.")

-          colnames(GRN@data$metadata)[which(colnames(GRN@data$metadata) == "sampleID")] = "sampleID_original"

-          

+        futile.logger::flog.warn("Renaming the first column to sampleID. However, this column already exists, it will be renamed accordingly.")

+        colnames(GRN@data$metadata)[which(colnames(GRN@data$metadata) == "sampleID")] = "sampleID_original"

+        

       } 

       colnames(GRN@data$metadata)[1] = "sampleID"

-

+      

       # Assume the ID is in column 1, has to be unique

       if (nrow(GRN@data$metadata) > length(unique(GRN@data$metadata$sampleID))) {

         message = paste0("The first column in the sample metadata table must contain only unique values, as it is used as sample ID. Make sure the values are unique.")

@@ -252,8 +252,8 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

       dplyr::mutate(has_RNA = sampleID  %in% samples_rna,

                     has_peaks = sampleID %in% samples_peaks,

                     has_both = has_RNA & has_peaks

-                    )

-

+      )

+    

     GRN@config$sharedSamples = dplyr::filter(GRN@data$metadata, has_both) %>% dplyr::pull(sampleID) %>% as.character()

     counts_peaks = as.data.frame(counts_peaks)

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

     # Store the raw peaks data efficiently as DESeq object only if it contains only integers, otherwise store as normal matrix

-

+    

     if (isIntegerMatrix(counts_peaks[, GRN@config$sharedSamples])) {

       GRN@data$peaks$counts_orig = DESeq2::DESeqDataSetFromMatrix(counts_peaks[, GRN@config$sharedSamples], 

-                                                           colData = dplyr::filter(GRN@data$metadata, has_both) %>% tibble::column_to_rownames("sampleID"),

-                                                           design = ~1)

+                                                                  colData = dplyr::filter(GRN@data$metadata, has_both) %>% tibble::column_to_rownames("sampleID"),

+                                                                  design = ~1)

     } else {

       GRN@data$peaks$counts_orig = as.matrix(counts_peaks[, GRN@config$sharedSamples])

@@ -276,14 +276,14 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     if (isIntegerMatrix(counts_rna[, GRN@config$sharedSamples])) {

       GRN@data$RNA$counts_orig = DESeq2::DESeqDataSetFromMatrix(counts_rna[, GRN@config$sharedSamples], 

-                                                          colData = dplyr::filter(GRN@data$metadata, has_both) %>% tibble::column_to_rownames("sampleID"),

-                                                          design = ~1)   

+                                                                colData = dplyr::filter(GRN@data$metadata, has_both) %>% tibble::column_to_rownames("sampleID"),

+                                                                design = ~1)   

     } else {

-        GRN@data$RNA$counts_orig = as.matrix(counts_rna[, GRN@config$sharedSamples])

+      GRN@data$RNA$counts_orig = as.matrix(counts_rna[, GRN@config$sharedSamples])

     }

-

-  

+    

+    

     # Consensus peaks

     GRN@data$peaks$consensusPeaks = .createConsensusPeaksDF(getCounts(GRN, type = "peaks", norm = TRUE, permuted = FALSE)) 

     GRN@data$peaks$consensusPeaks = GRN@data$peaks$consensusPeaks[match(getCounts(GRN, type = "peaks", norm = TRUE, permuted = FALSE)$peakID, GRN@data$peaks$consensusPeaks$peakID), ]

@@ -294,29 +294,29 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

     # Assume 0-based exclusive format, see https://arnaudceol.wordpress.com/2014/09/18/chromosome-coordinate-systems-0-based-1-based/ and http://genome.ucsc.edu/FAQ/FAQformat.html#format1 for details

     consensus.gr   = .constructGRanges(GRN@data$peaks$consensusPeaks, seqlengths = .getChrLengths(GRN@config$parameters$genomeAssembly), GRN@config$parameters$genomeAssembly, zeroBased = TRUE)

-

+    

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

     if (length(overlappingPeaks) > 0){

+      

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

+      

+      messageAll = paste0(" ", length(overlappingPeaks), 

+                          " overlapping peaks have been identified. The first ten are: ", paste0(ids[seq_len(min(10, length(ids)))], collapse = ","),

+                          ". This may not be what you want, since overlapping peaks may have a heigher weight in the network. "

+      )

+      

+      

+      if (allowOverlappingPeaks) {

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

-        

-        messageAll = paste0(" ", length(overlappingPeaks), 

-                            " overlapping peaks have been identified. The first ten are: ", paste0(ids[seq_len(min(10, length(ids)))], collapse = ","),

-                            ". This may not be what you want, since overlapping peaks may have a heigher weight in the network. "

-        )

-        

-        

-        if (allowOverlappingPeaks) {

-            

-            message = paste0(messageAll, "As allowOverlappingPeaks has been set to TRUE, this is only a warning and not an error.")

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

-        } else {

-            message = paste0(messageAll, "As allowOverlappingPeaks = FALSE (the default), this is an error and not a warning. You may want to regenerate the peak file, eliminate peak overlaps, and rerun this function")

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

-        }

-        

+        message = paste0(messageAll, "As allowOverlappingPeaks has been set to TRUE, this is only a warning and not an error.")

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

+      } else {

+        message = paste0(messageAll, "As allowOverlappingPeaks = FALSE (the default), this is an error and not a warning. You may want to regenerate the peak file, eliminate peak overlaps, and rerun this function")

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

+      }

+      

     }

   }

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

   # Add gene annotation once

   GRN = .populateGeneAnnotation(GRN)

-

+  

   GRN

 }

@@ -358,14 +358,14 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 }

 .removeScientificNotation_positions <- function(peakIDs.vec) {

-    ids = strsplit(peakIDs.vec, split = ":", fixed = TRUE)

-    ids_chr = sapply(ids, "[[", 1)

-    ids_pos = sapply(ids, "[[", 2)

-    ids_pos = strsplit(ids_pos, split = "-", fixed = TRUE)

-    start = sapply(ids_pos, "[[", 1)

-    end   = sapply(ids_pos, "[[", 2)

-    

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

+  ids = strsplit(peakIDs.vec, split = ":", fixed = TRUE)

+  ids_chr = sapply(ids, "[[", 1)

+  ids_pos = sapply(ids, "[[", 2)

+  ids_pos = strsplit(ids_pos, split = "-", fixed = TRUE)

+  start = sapply(ids_pos, "[[", 1)

+  end   = sapply(ids_pos, "[[", 2)

+  

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

 }

@@ -412,90 +412,90 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 }

 .populatePeakAnnotation <- function (GRN) {

+  

+  countsPeaks.clean = getCounts(GRN, type = "peaks", norm = TRUE, permuted = FALSE)

+  

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

+  

+  countsPeaks.m = as.matrix(dplyr::select(countsPeaks.clean, -peakID))

+  

+  rowMeans_peaks   = rowMeans(countsPeaks.m)

+  rowMedians_peaks = matrixStats::rowMedians(countsPeaks.m)

+  CV_peaks = matrixStats::rowSds(countsPeaks.m) /  rowMeans_peaks

+  

+  metadata_peaks = tibble::tibble(peak.ID = countsPeaks.clean$peakID, 

+                                  peak.mean = rowMeans_peaks, 

+                                  peak.median = rowMedians_peaks, 

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

+  } else {

-    countsPeaks.clean = getCounts(GRN, type = "peaks", norm = TRUE, permuted = FALSE)

-    

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

-    

-    countsPeaks.m = as.matrix(dplyr::select(countsPeaks.clean, -peakID))

-    

-    rowMeans_peaks   = rowMeans(countsPeaks.m)

-    rowMedians_peaks = matrixStats::rowMedians(countsPeaks.m)

-    CV_peaks = matrixStats::rowSds(countsPeaks.m) /  rowMeans_peaks

-    

-    metadata_peaks = tibble::tibble(peak.ID = countsPeaks.clean$peakID, 

-                                    peak.mean = rowMeans_peaks, 

-                                    peak.median = rowMedians_peaks, 

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

-    } else {

-        

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

-        genomeAssembly = GRN@config$parameters$genomeAssembly

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

-        consensusPeaks.gr  = .constructGRanges(consensusPeaks, seqlengths = .getChrLengths(genomeAssembly), genomeAssembly)

-        

-        # Add ChIPSeeker anotation

-        peaks.annotated = suppressMessages(ChIPseeker::annotatePeak(

-            consensusPeaks.gr,

-            tssRegion = c(-5000, 5000), # extended from 3kb to 5

-            TxDb = .getGenomeObject(genomeAssembly, type = "txbd"),

-            level = "gene", 

-            assignGenomicAnnotation = TRUE,  # the default

-            genomicAnnotationPriority = c("Promoter", "5UTR", "3UTR", "Exon", "Intron",

-                                          "Downstream", "Intergenic"),  # the default

-            annoDb = .getGenomeObject(genomeAssembly, type = "packageName"), # optional, if provided, extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added

-            sameStrand = FALSE, # the default

-            ignoreOverlap = FALSE, # the default

-            ignoreUpstream = FALSE, # the default

-            ignoreDownstream = FALSE, # the default

-            overlap = "TSS", # the default

-            verbose = TRUE # the default

-        ))

-        

-        GRN@annotation$consensusPeaks_obj = peaks.annotated

-        

-        peaks.annotated.df = as.data.frame(peaks.annotated)

-        peaks.annotated.df$annotation[grepl("Exon", peaks.annotated.df$annotation)] = "Exon"

-        peaks.annotated.df$annotation[grepl("Intron", peaks.annotated.df$annotation)] = "Intron"

-        

-        GRN@annotation$consensusPeaks = dplyr::left_join(GRN@annotation$consensusPeaks, 

-                                                         peaks.annotated.df  %>% 

-                                                             dplyr::select(peakID, annotation, tidyselect::starts_with("gene"), -geneId, distanceToTSS, ENSEMBL, SYMBOL, GENENAME) %>%

-                                                             dplyr::mutate(annotation  = as.factor(annotation), 

-                                                                           ENSEMBL = as.factor(ENSEMBL), 

-                                                                           GENENAME = as.factor(GENENAME),

-                                                                           SYMBOL = as.factor(SYMBOL)),

-                                                         by = c("peak.ID" = "peakID")) %>%

-            dplyr::rename(peak.gene.chr = geneChr,

-                          peak.gene.start = geneStart, 

-                          peak.gene.end = geneEnd, 

-                          peak.gene.length = geneLength, 

-                          peak.gene.strand = geneStrand, 

-                          peak.gene.name = GENENAME,

-                          peak.gene.distanceToTSS = distanceToTSS,

-                          peak.gene.ENSEMBL = ENSEMBL,

-                          peak.gene.symbol = SYMBOL,

-                          peak.annotation = annotation

-            )

-        

-        

-    }

-    

-    

-    

-    

-    # Also add GC content as annotation columns

-    GRN = .calcGCContentPeaks(GRN)

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

+    genomeAssembly = GRN@config$parameters$genomeAssembly

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

+    consensusPeaks.gr  = .constructGRanges(consensusPeaks, seqlengths = .getChrLengths(genomeAssembly), genomeAssembly)

+    

+    # Add ChIPSeeker anotation

+    peaks.annotated = suppressMessages(ChIPseeker::annotatePeak(

+      consensusPeaks.gr,

+      tssRegion = c(-5000, 5000), # extended from 3kb to 5

+      TxDb = .getGenomeObject(genomeAssembly, type = "txbd"),

+      level = "gene", 

+      assignGenomicAnnotation = TRUE,  # the default

+      genomicAnnotationPriority = c("Promoter", "5UTR", "3UTR", "Exon", "Intron",

+                                    "Downstream", "Intergenic"),  # the default

+      annoDb = .getGenomeObject(genomeAssembly, type = "packageName"), # optional, if provided, extra columns including SYMBOL, GENENAME, ENSEMBL/ENTREZID will be added

+      sameStrand = FALSE, # the default

+      ignoreOverlap = FALSE, # the default

+      ignoreUpstream = FALSE, # the default

+      ignoreDownstream = FALSE, # the default

+      overlap = "TSS", # the default

+      verbose = TRUE # the default

+    ))

+    

+    GRN@annotation$consensusPeaks_obj = peaks.annotated

+    

+    peaks.annotated.df = as.data.frame(peaks.annotated)

+    peaks.annotated.df$annotation[grepl("Exon", peaks.annotated.df$annotation)] = "Exon"

+    peaks.annotated.df$annotation[grepl("Intron", peaks.annotated.df$annotation)] = "Intron"

+    

+    GRN@annotation$consensusPeaks = dplyr::left_join(GRN@annotation$consensusPeaks, 

+                                                     peaks.annotated.df  %>% 

+                                                       dplyr::select(peakID, annotation, tidyselect::starts_with("gene"), -geneId, distanceToTSS, ENSEMBL, SYMBOL, GENENAME) %>%

+                                                       dplyr::mutate(annotation  = as.factor(annotation), 

+                                                                     ENSEMBL = as.factor(ENSEMBL), 

+                                                                     GENENAME = as.factor(GENENAME),

+                                                                     SYMBOL = as.factor(SYMBOL)),

+                                                     by = c("peak.ID" = "peakID")) %>%

+      dplyr::rename(peak.gene.chr = geneChr,

+                    peak.gene.start = geneStart, 

+                    peak.gene.end = geneEnd, 

+                    peak.gene.length = geneLength, 

+                    peak.gene.strand = geneStrand, 

+                    peak.gene.name = GENENAME,

+                    peak.gene.distanceToTSS = distanceToTSS,

+                    peak.gene.ENSEMBL = ENSEMBL,

+                    peak.gene.symbol = SYMBOL,

+                    peak.annotation = annotation

+      )

-    GRN

+  }

+  

+  

+  

+  

+  # Also add GC content as annotation columns

+  GRN = .calcGCContentPeaks(GRN)

+  

+  GRN

+  

 }

 .populateGeneAnnotation <- function (GRN) {

@@ -536,52 +536,52 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 #' @importFrom IRanges width

 #' @importFrom rlang .data `:=`

 .calcGCContentPeaks <- function(GRN) {

-    

-    futile.logger::flog.info(paste0("Calculate GC-content for peaks... "))

-    start = Sys.time()

-    genomeAssembly = GRN@config$parameters$genomeAssembly

-    #TODO: GC content for all peaks

-    genome = .getGenomeObject(genomeAssembly, type = "BSgenome")

-    

-    # Get peaks as GRanges object

-    query   = .constructGRanges(GRN@data$peaks$consensusPeaks, 

-                                seqlengths = .getChrLengths(genomeAssembly), 

-                                genomeAssembly)

-    

-    # Get DNAStringSet object

-    seqs_peaks = Biostrings::getSeq(genome, query)

-    

-    GC_content.df = Biostrings::letterFrequency(seqs_peaks, "GC") / Biostrings::letterFrequency(seqs_peaks, "ACGT")

-    

-    GC_content.df = GC_content.df %>%

-        tibble::as_tibble() %>%

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

-                      peak.ID = query$peakID,

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

-    

-    GC_classes.df = GC_content.df %>%

-        dplyr::group_by(GC_class) %>%

-        dplyr::summarise(n= dplyr::n(), peak_width_mean = mean(length), peak_width_sd = sd(length)) %>%

-        dplyr::ungroup() %>% 

-        tidyr::complete(GC_class, fill = list(n = 0)) %>%

-        dplyr::mutate(n_rel = .data$n / length(query))

-    

-    # TODO: Put where

-    #ggplot(GC_content.df, aes(GC.class)) + geom_histogram(stat = "count") + theme_bw()

-    

-    #ggplot(GC_classes.df , aes(GC.class, n_rel)) + geom_bar(stat = "identity") + theme_bw()

-    

-    GRN@annotation$consensusPeaks = dplyr::left_join(GRN@annotation$consensusPeaks, GC_content.df, by = "peak.ID") %>%

-        dplyr::rename( peak.GC.perc    = `G|C`,

-                       peak.width      = length,

-                       peak.GC.class   = GC_class)

-    

-    GRN@stats$peaks = list()

-    GRN@stats$peaks$GC = GC_classes.df

-    

-    .printExecutionTime(start)

-    

-    GRN

+  

+  futile.logger::flog.info(paste0("Calculate GC-content for peaks... "))

+  start = Sys.time()

+  genomeAssembly = GRN@config$parameters$genomeAssembly

+  #TODO: GC content for all peaks

+  genome = .getGenomeObject(genomeAssembly, type = "BSgenome")

+  

+  # Get peaks as GRanges object

+  query   = .constructGRanges(GRN@data$peaks$consensusPeaks, 

+                              seqlengths = .getChrLengths(genomeAssembly), 

+                              genomeAssembly)

+  

+  # Get DNAStringSet object

+  seqs_peaks = Biostrings::getSeq(genome, query)

+  

+  GC_content.df = Biostrings::letterFrequency(seqs_peaks, "GC") / Biostrings::letterFrequency(seqs_peaks, "ACGT")

+  

+  GC_content.df = GC_content.df %>%

+    tibble::as_tibble() %>%

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

+                  peak.ID = query$peakID,

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

+  

+  GC_classes.df = GC_content.df %>%

+    dplyr::group_by(GC_class) %>%

+    dplyr::summarise(n= dplyr::n(), peak_width_mean = mean(length), peak_width_sd = sd(length)) %>%

+    dplyr::ungroup() %>% 

+    tidyr::complete(GC_class, fill = list(n = 0)) %>%

+    dplyr::mutate(n_rel = .data$n / length(query))

+  

+  # TODO: Put where

+  #ggplot(GC_content.df, aes(GC.class)) + geom_histogram(stat = "count") + theme_bw()

+  

+  #ggplot(GC_classes.df , aes(GC.class, n_rel)) + geom_bar(stat = "identity") + theme_bw()

+  

+  GRN@annotation$consensusPeaks = dplyr::left_join(GRN@annotation$consensusPeaks, GC_content.df, by = "peak.ID") %>%

+    dplyr::rename( peak.GC.perc    = `G|C`,

+                   peak.width      = length,

+                   peak.GC.class   = GC_class)

+  

+  GRN@stats$peaks = list()

+  GRN@stats$peaks$GC = GC_classes.df

+  

+  .printExecutionTime(start)

+  

+  GRN

 }

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

@@ -637,8 +637,8 @@ filterData <- function (GRN,

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

       GRN@data$TFs$TF_peak_overlap[, "isFiltered"] = 0

     }

-  

-

+    

+    

     # Filter peaks

     futile.logger::flog.info("FILTER PEAKS")

     peakIDs.CV = .filterPeaksByMeanCV(GRN, 

@@ -660,7 +660,7 @@ filterData <- function (GRN,

     #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

-

+    

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

       GRN@data$TFs$TF_peak_overlap[, "isFiltered"] = as.integer (! rownames(GRN@data$TFs$TF_peak_overlap) %in% peaks_toKeep)

     }

@@ -671,9 +671,9 @@ filterData <- function (GRN,

     # Filter peaks

     futile.logger::flog.info("FILTER RNA-seq")

     genes.CV = .filterGenesByMeanCV(GRN, 

-                                      minMean = minNormalizedMeanRNA, maxMean = maxNormalizedMeanRNA, 

-                                      minCV = minCV_genes, maxCV = maxCV_genes) 

- 

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

     for (permutationCur in c(0)) {

       permIndex = as.character(permutationCur)

@@ -681,7 +681,7 @@ filterData <- function (GRN,

       GRN@data$RNA$counts_norm.l[[permIndex]]$isFiltered = rowMeans < minNormalizedMeanRNA 

     }

     nRowsFlagged = length(which(GRN@data$RNA$counts_norm.l[["0"]]$isFiltered))

-

+    

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

@@ -699,7 +699,7 @@ filterData <- function (GRN,

   startTime = Sys.time()

   if (is.null(minSize_peaks)) {

-      minSize_peaks = 1

+    minSize_peaks = 1

   }

   futile.logger::flog.info(paste0("Filter and sort peaks and remain only those on the following chromosomes: ", paste0(chrToKeep, collapse = ",")))

@@ -730,111 +730,111 @@ filterData <- function (GRN,

 }

 .filterPeaksByCV <- function (GRN, minCV = 0, maxCV = 1e7) {

+  

+  startTime = Sys.time()

+  

+  if (is.null(maxCV)) {

+    futile.logger::flog.info(paste0("Filter peaks by CV: Min CV = ", minCV))

-    startTime = Sys.time()

-    

-    if (is.null(maxCV)) {

-        futile.logger::flog.info(paste0("Filter peaks by CV: Min CV = ", minCV))

-        

-        peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, peak.CV >= minCV)

-        

-    } else {

-        futile.logger::flog.info(paste0("Filter peaks by CV: Min CV = ", minCV, ", max CV = ", maxCV))

-        

-        peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, peak.CV >= minCV, peak.CV <= maxCV)

-    }

-    

-    futile.logger::flog.info(paste0(" Number of peaks after filtering : ", nrow(peaksFiltered)))

+    peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, peak.CV >= minCV)

-    .printExecutionTime(startTime)

+  } else {

+    futile.logger::flog.info(paste0("Filter peaks by CV: Min CV = ", minCV, ", max CV = ", maxCV))

-    peaksFiltered$peak.ID

+    peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, peak.CV >= minCV, peak.CV <= maxCV)

+  }

+  

+  futile.logger::flog.info(paste0(" Number of peaks after filtering : ", nrow(peaksFiltered)))

+  

+  .printExecutionTime(startTime)

+  

+  peaksFiltered$peak.ID

 }

 .filterPeaksByMeanCV <- function (GRN, minMean = 0, maxMean = NULL, minCV = 0, maxCV = NULL) {

+  

+  startTime = Sys.time()

+  

+  futile.logger::flog.info(paste0(" Number of peaks before filtering : ", nrow(GRN@annotation$consensusPeaks)))

+  

+  if (is.null(minCV)) {

+    minCV = 0

+  }

+  

+  if (is.null(maxCV)) {

+    futile.logger::flog.info(paste0("  Filter peaks by CV: Min = ", minCV))

+    maxCV = 9e+99

-    startTime = Sys.time()

-    

-    futile.logger::flog.info(paste0(" Number of peaks before filtering : ", nrow(GRN@annotation$consensusPeaks)))

-    

-    if (is.null(minCV)) {

-        minCV = 0

-    }

-    

-    if (is.null(maxCV)) {

-        futile.logger::flog.info(paste0("  Filter peaks by CV: Min = ", minCV))

-        maxCV = 9e+99

-        

-    } else {

-        futile.logger::flog.info(paste0("  Filter peaks by CV: Min = ", minCV, ", Max = ", maxCV))

-    }

-    

-    

-    if (is.null(minMean)) {

-        minMean = 0

-    }

-    

-    if (is.null(maxMean)) {

-        futile.logger::flog.info(paste0("  Filter peaks by mean: Min = ", minMean))

-        maxMean = 9e+99

-    } else {

-        futile.logger::flog.info(paste0("  Filter peaks by mean: Min = ", minMean, ", Max = ", maxMean))  

-    }   

-    

-    

-    peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, 

-                                  peak.CV >= minCV, peak.CV <= maxCV, 

-                                  peak.mean >= minMean, peak.mean <= maxMean)

-    

-    futile.logger::flog.info(paste0(" Number of peaks after filtering : ", nrow(peaksFiltered)))

-    

-    .printExecutionTime(startTime)

-    

-    peaksFiltered$peak.ID

+  } else {

+    futile.logger::flog.info(paste0("  Filter peaks by CV: Min = ", minCV, ", Max = ", maxCV))

+  }

+  

+  

+  if (is.null(minMean)) {

+    minMean = 0

+  }

+  

+  if (is.null(maxMean)) {

+    futile.logger::flog.info(paste0("  Filter peaks by mean: Min = ", minMean))

+    maxMean = 9e+99

+  } else {

+    futile.logger::flog.info(paste0("  Filter peaks by mean: Min = ", minMean, ", Max = ", maxMean))  

+  }   

+  

+  

+  peaksFiltered = dplyr::filter(GRN@annotation$consensusPeaks, 

+                                peak.CV >= minCV, peak.CV <= maxCV, 

+                                peak.mean >= minMean, peak.mean <= maxMean)

+  

+  futile.logger::flog.info(paste0(" Number of peaks after filtering : ", nrow(peaksFiltered)))

+  

+  .printExecutionTime(startTime)

+  

+  peaksFiltered$peak.ID

 }

 .filterGenesByMeanCV <- function (GRN, minMean = 0, maxMean = NULL, minCV = 0, maxCV = NULL) {

+  

+  startTime = Sys.time()

+  

+  futile.logger::flog.info(paste0(" Number of genes before filtering : ", nrow(GRN@annotation$genes)))

+  

+  if (is.null(minCV)) {

+    minCV = 0

+  }

+  

+  if (is.null(maxCV)) {

+    futile.logger::flog.info(paste0("  Filter genes by CV: Min = ", minCV))

+    maxCV = 9e+99

-    startTime = Sys.time()

-    

-    futile.logger::flog.info(paste0(" Number of genes before filtering : ", nrow(GRN@annotation$genes)))

-    

-    if (is.null(minCV)) {

-        minCV = 0

-    }

-    

-    if (is.null(maxCV)) {

-        futile.logger::flog.info(paste0("  Filter genes by CV: Min = ", minCV))

-        maxCV = 9e+99

-        

-    } else {

-        futile.logger::flog.info(paste0("  Filter genes by CV: Min = ", minCV, ", Max = ", maxCV))

-    }

-    

-    

-    if (is.null(minMean)) {

-        minMean = 0

-    }

-    

-    

-    if (is.null(maxMean)) {

-        futile.logger::flog.info(paste0("  Filter genes by mean: Min = ", minMean))

-        maxMean = 9e+99

-    } else {

-        futile.logger::flog.info(paste0("  Filter genes by mean: Min = ", minMean, ", Max = ", maxMean))  

-    }   

-    

-    

-    genesFiltered = dplyr::filter(GRN@annotation$genes, 

-                                  gene.CV >= minCV, gene.CV <= maxCV, 

-                                  gene.mean >= minMean, gene.mean <= maxMean)

-

-    

-    futile.logger::flog.info(paste0(" Number of genes after filtering : ", nrow(genesFiltered)))

-    

-    .printExecutionTime(startTime)

-    

-    genesFiltered$gene.ENSEMBL 

+  } else {

+    futile.logger::flog.info(paste0("  Filter genes by CV: Min = ", minCV, ", Max = ", maxCV))

+  }

+  

+  

+  if (is.null(minMean)) {

+    minMean = 0

+  }

+  

+  

+  if (is.null(maxMean)) {

+    futile.logger::flog.info(paste0("  Filter genes by mean: Min = ", minMean))

+    maxMean = 9e+99

+  } else {

+    futile.logger::flog.info(paste0("  Filter genes by mean: Min = ", minMean, ", Max = ", maxMean))  

+  }   

+  

+  

+  genesFiltered = dplyr::filter(GRN@annotation$genes, 

+                                gene.CV >= minCV, gene.CV <= maxCV, 

+                                gene.mean >= minMean, gene.mean <= maxMean)

+  

+  

+  futile.logger::flog.info(paste0(" Number of genes after filtering : ", nrow(genesFiltered)))

+  

+  .printExecutionTime(startTime)

+  

+  genesFiltered$gene.ENSEMBL 

 }

@@ -855,9 +855,9 @@ filterData <- function (GRN,

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

-    

+  

   checkmate::assertClass(GRN, "GRN")

   checkmate::assertDirectoryExists(motifFolder)

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

@@ -867,12 +867,12 @@ addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPatte

   checkmate::assertFlag(forceRerun)

   if (is.null(GRN@data$TFs$translationTable) | is.null(GRN@data$TFs$translationTable) | is.null(GRN@config$allTF)  | is.null(GRN@config$directories$motifFolder) | forceRerun) {

-

+    

     GRN@config$TFBS_fileEnding  = fileEnding

     GRN@config$TFBS_filePattern = filesTFBSPattern

     GRN@data$TFs$translationTable = .getFinalListOfTFs(motifFolder, filesTFBSPattern, fileEnding, TFs, nTFMax, getCounts(GRN, type = "rna", norm = TRUE, permuted = FALSE))

-

+    

     GRN@data$TFs$translationTable = GRN@data$TFs$translationTable %>%

       dplyr::select(c("HOCOID", "ENSEMBL")) %>%

       dplyr::mutate(TF.name = HOCOID, TF.ENSEMBL = ENSEMBL) 

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

     # GRN@config$TF_list = list()

     # GRN@config$TF_list[["all_TFBS"]] =GRN@config$allTF

     GRN@config$directories$motifFolder = motifFolder

-  

-  

+    

+    

   } 

   GRN

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

   files = .createFileList(folder_input_TFBS, "*.bed*", recursive = FALSE, ignoreCase = FALSE, verbose = FALSE)

   TFsWithTFBSPredictions = gsub(pattern = filesTFBSPattern, "", tools::file_path_sans_ext(basename(files), compression = TRUE))

   TFsWithTFBSPredictions = gsub(pattern = fileEnding, "", TFsWithTFBSPredictions)

-

+  

   futile.logger::flog.info(paste0("Found ", length(TFsWithTFBSPredictions), " matching TFs: ", paste0(TFsWithTFBSPredictions, collapse = ", ")))

-

+  

   # Filter TFs

   if (length(TFs) == 1 && TFs == "all") {

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

       allTF = allTF[seq_len(nTFMax)]

       futile.logger::flog.info(paste0("Updated list of TFs: ", paste0(allTF, collapse = ", ")))

     } 

-

+    

   }

   nTF = length(allTF)

@@ -974,18 +974,18 @@ addTFBS <- function(GRN, motifFolder, TFs = "all", nTFMax = NULL, filesTFBSPatte

 #' GRN = overlapPeaksAndTFBS(GRN, nCores = 2, forceRerun = FALSE)

 #' @export

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

-

+  

   GRN = .addFunctionLogToObject(GRN)

-    

+  

   checkmate::assertClass(GRN, "GRN")

   checkmate::assertIntegerish(nCores, lower = 1)

   checkmate::assertFlag(forceRerun)

   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..."))

-

+