Various changes, restructuring

src/GRaNIE/DESCRIPTION

 Package: GRaNIE
 Title: GRaNIE: Reconstruction cell type specific gene regulatory networks including enhancers using chromatin accessibility and RNA-seq data
-Version: 0.14.6
+Version: 0.17.0
 Encoding: UTF-8
 Authors@R: c(person("Christian", "Arnold", email =
-        "christian.arnold@embl.de", role = c("cre","aut")),
+        "chrarnold@web.de", role = c("cre","aut")),
         person("Judith", "Zaugg", email =
         "judith.zaugg@embl.de", role = c("aut")),
         person("Rim", "Moussa", email =
@@ -63,7 +63,7 @@ Imports:
     magrittr,
     tibble
 Depends:
-    R (>= 3.6),
+    R (>= 4.1.0),
     tidyverse,
     topGO
 Suggests:
@@ -77,6 +77,10 @@ Suggests:
     TxDb.Mmusculus.UCSC.mm10.knownGene,
     TxDb.Mmusculus.UCSC.mm9.knownGene,
     IHW,
+    biomaRt,
+    clusterProfiler,
+    ReactomePA,
+    DOSE,
     ChIPseeker,
     testthat (>= 3.0.0)
 VignetteBuilder: knitr
@@ -90,9 +94,12 @@ biocViews:
   Genetics,
   Transcriptomics,
   ATACSeq,
-  RNASeq
+  RNASeq,
+  Network,
+  Transcription,
+  ChIPSeq
 License: LGPL (>= 3)
-LazyData: true
+LazyData: false
 URL: https://grp-zaugg.embl-community.io/GRaNIE
 BugReports: https://git.embl.de/grp-zaugg/GRaNIE/issues
 RoxygenNote: 7.1.2

src/GRaNIE/NAMESPACE

@@ -7,6 +7,7 @@ export(addData)
 export(addData_TFActivity)
 export(addTFBS)
 export(add_TF_gene_correlation)
+export(build_eGRN_graph)
 export(calculateCommunitiesEnrichment)
 export(calculateCommunitiesStats)
 export(calculateGeneralEnrichment)
@@ -34,12 +35,17 @@ export(plotGeneralGraphStats)
 export(plotPCA_all)
 export(plotTFEnrichment)
 export(plot_stats_connectionSummary)
+import(BiocManager)
+import(DOSE)
 import(GenomicRanges)
 import(IRanges)
+import(ReactomePA)
 import(checkmate)
+import(clusterProfiler)
 import(ggplot2)
 import(graphics)
 import(grid)
+import(magrittr)
 import(patchwork)
 import(rlang)
 import(tibble)

src/GRaNIE/R/DataClasses.R

@@ -44,7 +44,7 @@
 #' @slot stats Stores statistical and summary information for a \code{GRN} network. Currently, connection details are stored here.
 #' 
 #' @slot visualization Stores visualization results, currently always empty. Feature in development.
-#' @slot isDev Flag whether this is an object of class \code{GRNdev}.
+#' @slot isDev Flag whether this is an object from the development version of the package
 
 #' @keywords GRN-class, GRN
 #' @name GRN-class
@@ -57,6 +57,7 @@ setClass("GRN", representation = representation(
   data          = "list",
   stats         = "list",
   visualization = "list",
+  graph         = "list",
   isDev         = "logical"
   )
 )
@@ -521,7 +522,7 @@ setMethod("show",
             if (!is.null(GRN@data$peaks$consensusPeaks)) {
               nPeaks_filt  = nPeaks(GRN, filter = TRUE)
               nPeaks_all   = nPeaks(GRN, filter = FALSE)
-              cat(" Number of peaks (filtered, all): ", nPeaks_filt, ", ",  nPeaks_all, "\n")
+              cat(" Number of peaks (filtered, all): ", nPeaks_filt, ", ",  nPeaks_all, "\n", sep = "")
             } else {
               cat (" Number of peaks: No peak data found.\n")
             }
@@ -537,7 +538,7 @@ setMethod("show",
                 
               }
              
-              cat(" Number of genes (filtered, all): ", nGenes_filt, ", ",  nGenes_all, "\n")
+              cat(" Number of genes (filtered, all): ", nGenes_filt, ", ",  nGenes_all, "\n", sep = "")
             } else {
               cat (" Number of genes: No RNA-seq data found.\n")
             }
@@ -557,21 +558,84 @@ setMethod("show",
             
             cat("Connections:\n")
             if (!is.null(GRN@connections)) {
-              
-              nGenes_filt = nGenes(GRN, filter = TRUE)
-              nGenes_all  = nGenes(GRN, filter = FALSE)
-              
-              if (checkmate::testClass(GRN@data$RNA$counts_orig, "DESeqDataSet")) {
+
+                if (!is.null(GRN@connections$TF_peaks$`0`$main)) {
+                    maxFDR = round(max(GRN@connections$TF_peaks$`0`$main$TF_peak.fdr), 1)
+                    cat(" TF-peak links (", nrow(GRN@connections$TF_peaks$`0`$main), " with FDR < ", maxFDR, ")\n", sep = "")
+                } else {
+                    cat(" TF-peak links: none found\n")
+                }
+                if (!is.null(GRN@connections$peak_genes$`0`)) {
+                    cat(" peak-gene links (", nrow(GRN@connections$peak_genes$`0`)," with ", GRN@config$parameters$promoterRange, " bp promoter range)\n", sep = "")
+                } else {
+                    cat(" peak-gene links: none found\n")
+                }
                 
-              } else {
                 
-              }
+                if (!is.null(GRN@connections$all.filtered$`0`)) {
+                    nRows = nrow(GRN@connections$all.filtered$`0`)
+                    max_TF_peak_FDR = round(max(GRN@connections$all.filtered$`0`$TF_peak.fdr), 1)
+                    max_peak_gene_FDR = round(max(GRN@connections$all.filtered$`0`$peak_gene.p_adj), 1)
+                    cat(" TF-peak-gene links (", nRows," with TF-peak FDR ", max_TF_peak_FDR, " and peak-gene FDR ",  max_peak_gene_FDR , ")\n", sep = "")
+                    
+                } else {
+                    cat(" TF-peak-gene links (filtered): none found\n")
+                }
+            
+
               
-              cat(" Number of genes (filtered, all): ", nGenes_filt, ", ",  nGenes_all, "\n")
             } else {
-              cat (" Number of genes: No RNA-seq data found.\n")
+              cat (" Connections: none found")
             }
             
+            cat("Network-related:\n")
+            
+            if (!is.null(GRN@graph$TF_gene)) {
+                
+                cat(" eGRN network:\n")
+                nNodes = length(igraph::V(GRN@graph$TF_gene$graph))
+                nEdges = length(igraph::E(GRN@graph$TF_gene$graph))
+                cat("   TF-gene eGRN: ", nNodes, " nodes and ", nEdges, " edges\n",sep = "")
+ 
+                nNodes = length(igraph::V(GRN@graph$TF_peak_gene$graph))
+                nEdges = length(igraph::E(GRN@graph$TF_peak_gene$graph))
+                cat("   TF-peak-gene eGRN: ", nNodes, " nodes and ", nEdges, " edges\n",sep = "")              
+                
+                # Community identification (no, yes and how many and how many nodes each)
+                cat(" Communities (TF-gene):\n")
+                df = igraph::vertex.attributes(GRN@graph[["TF_gene"]]$graph)
+                if (!is.null(df)) {
+                    communities = df %>% as.data.frame() %>% dplyr::count(community) %>% dplyr::arrange(desc(n))
+                    cat("  Communities, sorted by size (n = Number of nodes): ", paste0(communities$community, " (n=", communities$n, collapse = "), "), ")\n", sep = "")
+                } else {
+                    cat("  None found\n")
+                }
+                
+                cat(" Enrichment (TF-gene):\n")
+                if (!is.null(GRN@stats$Enrichment$byTF)) {
+                    ontologies = names(GRN@stats$Enrichment$general)
+                    cat("  Overall network: ", paste0(ontologies, collapse = " & "), "\n",sep = "")
+                } else {
+                    cat("  Overall network: no enrichment results found\n")
+                }
+                
+                if (!is.null(GRN@stats$Enrichment$byTF)) {
+                    ontologies = names(GRN@stats$Enrichment$byCommunity[[1]])
+                    cat("  Communities: ", paste0(ontologies, collapse = " & "), "\n",sep = "")
+                } else {
+                    cat("  Communities: no enrichment results found\n")
+                }
+                
+                if (!is.null(GRN@stats$Enrichment$byTF)) {
+                    ontologies = names(GRN@stats$Enrichment$byTF[[1]])
+                    cat("  TFs: ", paste0(ontologies, collapse = " & "), "\n",sep = "")
+                } else {
+                    cat("  TFs: no enrichment results found\n")
+                }
+                
+            } else {
+                cat("  eGRN network: not found\n")
+            }
             
           }
 )

src/GRaNIE/R/diffTF_sharedFunctions.R

@@ -41,8 +41,7 @@
 .normalizeCounts <- function(rawCounts, method = "quantile", idColumn, removeCols = c(), returnDESeqObj = FALSE) {
   
   start = Sys.time()
-  .checkAndLoadPackages(c("tidyverse", "futile.logger", "checkmate", "limma"), verbose = FALSE)
-  
+
   futile.logger::flog.info(paste0("Normalize counts. Method: ", method, ", ID column: ", idColumn))
   
   if (is.null(idColumn)) {
@@ -276,8 +275,6 @@
 
 .finalizeClassificationAndAppend <- function(output.global.TFs, median.cor.tfs, act.rep.thres.l, par.l, t.cor.sel.matrix, t.cor.sel.matrix.non, significanceThreshold_Wilcoxon = 0.05) {
     
-   .checkAndLoadPackages(c("progress"), verbose = FALSE)
-  
     start = Sys.time()
     futile.logger::flog.info(paste0("Finalize classification"))
     colnameMedianCor         = paste0("median.cor.tfs")

src/GRaNIE/R/helperFunctions.R

-
-.checkAndLoadPackages <- function(packages, verbose = FALSE) {
-  
-  .checkAndInstallMissingPackages(packages, verbose = verbose)
-    
-    if ("tidyverse" %in% packages | "readr" %in% packages) {
-        if (packageVersion("readr") < "1.2.0") {
-            message = paste0("Version of readr library is too old, at least version 1.2.0 is required). Execute the following in R: install.packages('readr') ") 
-            .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)
-        }
-    }
-  
-  for (packageCur in packages) {
-    suppressMessages(library(packageCur, character.only = TRUE))
-  }
-  
-}
-
 .printMultipleGraphsPerPage <- function(plots.l, nCol = 1, nRow = 1, pdfFile = NULL, height = NULL, width = NULL, verbose = FALSE) {
     
-    #.checkAndLoadPackages(c("grDevices", "gridExtra"), verbose = verbose)  
-    
     checkmate::assertList(plots.l, min.len = 1)
     checkmate::assertInt(nCol, lower = 1)
     checkmate::assertInt(nRow, lower = 1)
@@ -83,8 +63,6 @@
 
 .startLogger <- function(logfile, level, removeOldLog = TRUE, appenderName = "consoleAndFile", verbose = TRUE) {
   
-  #.checkAndLoadPackages(c("futile.logger"), verbose = FALSE)
-  
   checkmate::assertSubset(level, c("TRACE", "DEBUG", "INFO", "WARN", "ERROR", "FATAL"))
   checkmate::assertFlag(removeOldLog)
   checkmate::assertSubset(appenderName, c("console", "file", "consoleAndFile"))
@@ -114,7 +92,6 @@
 
 .printParametersLog <- function(par.l, verbose = FALSE) {
   
-  #.checkAndLoadPackages(c("futile.logger"), verbose = verbose)  
   checkmate::assertList(par.l)
   futile.logger::flog.info(paste0("PARAMETERS:"))
   for (parCur in names(par.l)) {
@@ -128,55 +105,11 @@
 # PACKAGE LOADING AND DETACHING FUNCTIONS #
 ###########################################
 
-.checkAndInstallMissingPackages <- function(packages.vec, verbose = FALSE) {
-  
-  
-  if (verbose) cat(paste0("Trying to automatically install missing packages. If this fails, install them manually..."))
-  
-  packagesToInstall = setdiff(packages.vec, rownames(installed.packages()))
-  
-  
-  if (length(packagesToInstall) > 0) {
-    if (verbose) cat(paste0("Could not find the following packages: ", paste( packagesToInstall , collapse = ", ")))
-    install.packages(packagesToInstall, repos = "http://cran.rstudio.com/")  
-    
-    if (getRversion() < "3.5.0") { 
-        
-        source("http://bioconductor.org/biocLite.R")
-        for (packageCur in packagesToInstall) {
-            biocLite(packageCur, suppressUpdates = TRUE, ask = FALSE)
-        }
-    } else {
-        
-        if (!requireNamespace("BiocManager", quietly = TRUE))
-            install.packages("BiocManager")
-        
-        for (packageCur in packagesToInstall) {
-            BiocManager::install(packageCur, update = FALSE, ask = FALSE)
-        }
+.checkOntologyPackageInstallation <- function(pkg) {
+    if (!is.installed(pkg)) {
+        message = paste0("The package ", pkg, " is not installed, which is however needed for the chosen ontology enrichment. Please install it and re-run this function or change the ontology.")
+        .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)
     }
-    
-   
-    
-    
-    
-  } else {
-    if (verbose) cat("All packages are already installed")
-  }
-  
-
-}
-
-.detachAllPackages <- function() {
-  
-  basic.packages <- c("package:stats","package:graphics","package:grDevices","package:utils","package:datasets","package:methods","package:base")
-  
-  package.list <- search()[ifelse(unlist(gregexpr("package:",search())) == 1,TRUE,FALSE)]
-  
-  package.list <- setdiff(package.list,basic.packages)
-  
-  if (length(package.list) > 0)  for (package in package.list) detach(package, character.only = TRUE)
-  
 }
 
 .clearOpenDevices <- function() {
@@ -254,7 +187,6 @@
 
 .initBiocParallel <- function(nWorkers, verbose = FALSE) {
   
-  #.checkAndLoadPackages(c("BiocParallel"), verbose = verbose)  
   checkmate::assertInt(nWorkers, lower = 1)    
   checkmate::assertInt(BiocParallel::multicoreWorkers())
   
@@ -270,7 +202,6 @@
 
 .execInParallelGen <- function(nCores, returnAsList = TRUE, listNames = NULL, iteration, abortIfErrorParallel = TRUE, verbose = TRUE, functionName, ...) {
   
-  #.checkAndLoadPackages(c("BiocParallel"), verbose = verbose)  
   start.time  <-  Sys.time()
   
   checkmate::assertInt(nCores, lower = 1)
@@ -414,7 +345,7 @@
   gr = GenomicRanges::makeGRangesFromDataFrame(annotation, keep.extra.columns = TRUE, seqinfo = seqlengths, starts.in.df.are.0based = zeroBased, ...)
   
   # Check whether there are out-of-bound sequences, and abort if there are. This should not happen
-  seq_outOfBound = which(end(gr) > seqlengths(gr)[as.character(GenomeInfoDb::seqnames(gr))])
+  seq_outOfBound = which(end(gr) > GenomeInfoDb::seqlengths(gr)[as.character(GenomeInfoDb::seqnames(gr))])
   if (length(seq_outOfBound) > 0) {
     
       index = seq_len(min(length(seq_outOfBound),10))
@@ -447,6 +378,17 @@
   gr
 }
 
+.checkAndLoadPackages <- function(packageList, verbose = FALSE) {
+    
+    for (packageCur in packageList) {
+        if (!is.installed(packageCur)) {
+            message = paste0("The suggested package \"", packageCur, "\" is not yet installed, but it is required for this function and in this context. Please install it and re-run this function.")
+            .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)
+        }
+    }
+    
+}
+
 .getGenomeObject <- function(genomeAssembly, type = "txbd") {
     
     
@@ -515,10 +457,8 @@
 
 .checkAndLoadPackagesGenomeAssembly <- function(genomeAssembly) {
     
-    #.checkAndLoadPackages("checkmate")
     checkmate::assertSubset(genomeAssembly, c("hg19","hg38", "mm9", "mm10"))
     
-    
     if (genomeAssembly == "hg38") {
         
         packages = c("org.Hs.eg.db", "TxDb.Hsapiens.UCSC.hg38.knownGene")
@@ -609,17 +549,6 @@
     
 }
 
-loadAllPackages <- function() {
-  
-  # .initFunctionsScript(packagesReq = NULL, minRVersion = "3.5.0", warningsLevel = 1, disableScientificNotation = TRUE)
-  .checkAndLoadPackages(c("tidyverse", "futile.logger", "checkmate", "reshape2", "data.table", "rlist", "matrixStats", "Matrix",
-                         "bedr", "GenomicRanges", 
-                         "RColorBrewer", "igraph", "colorspace", "pheatmap",
-                         "preprocessCore", "DESeq2", "csaw", "robust", "PCAtools", "progress"), verbose = FALSE
-  )
-  
-}
-
 
 .asSparseMatrix <- function(matrix, convertNA_to_zero=TRUE, dimnames = NULL) {
   

src/GRaNIE/docs/404.html

@@ -49,7 +49,7 @@
       </button>
       <span class="navbar-brand">
         <a class="navbar-link" href="https://grp-zaugg.embl-community.io/GRaNIE/index.html">GRaNIE</a>
-        <span class="version label label-danger" data-toggle="tooltip" data-placement="bottom" title="Unreleased version">0.14.6</span>
+        <span class="version label label-danger" data-toggle="tooltip" data-placement="bottom" title="Unreleased version">0.17.0</span>
       </span>
     </div>