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

+Version: 0.99.4

 Encoding: UTF-8

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

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

@@ -370,23 +370,51 @@ addData <- function(GRN, counts_peaks, normalization_peaks = "DESeq_sizeFactor",

 }

-.loadPrecompiledAnnotationData <- function(genomeAssembly) {

-  

-  futile.logger::flog.info(paste0("Reading pre-compiled genome annotation data "))

-  

-  if (is.null(geneAnnotation[[genomeAssembly]])) {

+.retrieveAnnotationData <- function(genomeAssembly) {

-    message = "The genome version you specified is not contained in the pre-compiled genome annotation list. Currently, only hg19, hg38 and mm10 are supported. Contact us."

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

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

+    

+    params.l = .getBiomartParameters(genomeAssembly)

+    

+    columnsToRetrieve = c("chromosome_name", "start_position", "end_position",

+                          "strand", "ensembl_gene_id", "gene_biotype", "external_gene_name")

+    

+    geneAnnotation <- NULL

+    attempt <- 0

+    mirrors = c('www', 'uswest', 'useast', 'asia')

+    while(!is.data.frame(geneAnnotation) && attempt <= 3 ) {

+        attempt <- attempt + 1

+        geneAnnotation = tryCatch({ 

+            ensembl = biomaRt::useEnsembl(biomart = "genes", host = params.l[["host"]], dataset = params.l[["dataset"]], mirror = mirrors[attempt])

+            biomaRt::getBM(attributes = columnsToRetrieve, mart = ensembl)

+            

+        }

+        )

+    } 

+    

+    

+    if (!is.data.frame(geneAnnotation)) {

+        

+        error_Biomart = "A temporary error occured with biomaRt::getBM or biomaRt::useEnsembl. This is often caused by an unresponsive Ensembl site. Try again at a later time. Note that this error is not caused by GRaNIE but external services."

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

+        return(NULL)

+        

+    }

+    

+    geneAnnotation %>%

+        as_tibble() %>%

+        dplyr::filter(stringr::str_length(chromosome_name) <= 5) %>%

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

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

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

-  } 

-  

-  geneAnnotation[[genomeAssembly]]

 }

 .getAllGeneTypesAndFrequencies <- function(genomeAssembly, verbose = TRUE) {

-  geneAnnotation.df = .loadPrecompiledAnnotationData(genomeAssembly)

+  geneAnnotation.df = .retrieveAnnotationData(genomeAssembly)

   return(table(geneAnnotation.df$gene.type))

 }

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

     stop("Not yet implemented")

   }

-  genes.filt.df = .loadPrecompiledAnnotationData(GRN@config$parameters$genomeAssembly)

+  genes.filt.df = .retrieveAnnotationData(GRN@config$parameters$genomeAssembly)

   if (!is.null(gene.types)) {

     if (! "all" %in% gene.types) {

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

   rowMedians_rna = matrixStats::rowMedians(countsRNA.m)

   CV_rna = matrixStats::rowSds(countsRNA.m) /  rowMeans_rna

-  # This object is internally loaded and lives inside the R directory within sysdata.rda

-  # load("R/sysdata.rda")

-  genomeAnnotation.df = geneAnnotation[[GRN@config$parameters$genomeAssembly]]

+  genomeAnnotation.df = .retrieveAnnotationData(GRN@config$parameters$genomeAssembly)

   metadata_rna = tibble::tibble(gene.ENSEMBL = countsRNA.clean$ENSEMBL, 

                                 gene.mean = rowMeans_rna, 

@@ -2431,11 +2457,11 @@ addConnections_peak_gene <- function(GRN, overlapTypeGene = "TSS", corMethod = "

     end(subject.gr) = start(subject.gr)

   } 

-  overlapsAll = GenomicRanges::findOverlaps(query, subject.gr, 

+  overlapsAll = suppressWarnings(GenomicRanges::findOverlaps(query, subject.gr, 

                                             minoverlap=1,

                                             type="any",

                                             select="all",

-                                            ignore.strand=TRUE)

+                                            ignore.strand=TRUE))

   query_row_ids  = S4Vectors::queryHits(overlapsAll)

   subject_rowids = S4Vectors::subjectHits(overlapsAll)

@@ -2596,7 +2622,7 @@ addConnections_peak_gene <- function(GRN, overlapTypeGene = "TSS", corMethod = "

   # if a peak overlaps with a gene, should the same gene be reported as the connection?

   # OVERLAP OF PEAKS AND EXTENDED GENES

-  overlaps.sub.df = .calculatePeakGeneOverlaps(GRN, consensusPeaks, peak.TADs.df, neighborhoodSize = neighborhoodSize, 

+  overlaps.sub.df = .calculatePeakGeneOverlaps(GRN, allPeaks = consensusPeaks, peak.TADs.df, neighborhoodSize = neighborhoodSize, 

                                                genomeAssembly = genomeAssembly, gene.types = gene.types, overlapTypeGene = overlapTypeGene) 

@@ -4388,3 +4414,110 @@ changeOutputDirectory <- function(GRN, outputDirectory = ".") {

 }

+

+.createTables_peakGeneQC <- function(peakGeneCorrelations.all.cur, networkType_details, colors_vec, range) {

+    

+    d = peakGeneCorrelations.all.cur %>% 

+        dplyr::group_by(r_positive, class, peak_gene.p.raw.class) %>%

+        dplyr::summarise(n = dplyr::n()) %>%

+        dplyr::ungroup()

+    

+    # Some classes might be missing, add them here with explicit zeros

+    for (r_pos in c(TRUE, FALSE)) {

+        for (classCur in networkType_details){

+            for (pclassCur in levels(peakGeneCorrelations.all.cur$peak_gene.p.raw.class)) {

+                

+                row = which(d$r_positive == r_pos & d$class == classCur & as.character(d$peak_gene.p.raw.class) == as.character(pclassCur))

+                if (length(row) == 0) {

+                    d = tibble::add_row(d, r_positive = r_pos, class = classCur, peak_gene.p.raw.class = pclassCur, n = 0)

+                }

+            }

+        }

+    }

+    

+    # Restore the "ordered" factor for class

+    d$peak_gene.p.raw.class = factor(d$peak_gene.p.raw.class, ordered = TRUE, levels =  levels(peakGeneCorrelations.all.cur$peak_gene.p.raw.class))

+    

+    

+    # Normalization factors

+    dsum = d %>%

+        dplyr::group_by(r_positive, class) %>%

+        dplyr::summarise(sum_n = sum(.data$n))

+    

+    

+    # Summarize per bin

+    d2 = d %>%

+        dplyr::group_by(class, peak_gene.p.raw.class)%>%

+        dplyr::summarise(sum_pos = .data$n[r_positive],

+                         sum_neg = .data$n[!r_positive]) %>%

+        dplyr::mutate(ratio_pos_raw = sum_pos / sum_neg,

+                      enrichment_pos = sum_pos / (sum_pos + sum_neg),

+                      ratio_neg = 1 - enrichment_pos) %>%

+        dplyr::ungroup()

+    

+    # Compare between real and permuted

+    normFactor_real = dplyr::filter(dsum, class ==  !! (networkType_details[1])) %>%  dplyr::pull(sum_n) %>% sum() /

+        dplyr::filter(dsum, class ==  !! (networkType_details[2])) %>%  dplyr::pull(sum_n) %>% sum()

+    

+    # ratio_norm not used currently, no normalization necessary here or not even useful because we dont want to normalize the r_pos and r_neg ratios: These are signal in a way. Only when comparing between real and permuted, we have to account for sample size for corrections

+    d3 = d %>%

+        dplyr::group_by(peak_gene.p.raw.class, r_positive) %>%

+        dplyr::summarise(n_real     = .data$n[class == !! (names(colors_vec)[1]) ],

+                         n_permuted = .data$n[class == !! (names(colors_vec)[2]) ]) %>%

+        dplyr::ungroup() %>%

+        dplyr::mutate(ratio_real_raw = n_real / n_permuted,

+                      ratio_real_norm = ratio_real_raw / normFactor_real,

+                      enrichment_real      = n_real / (n_real + n_permuted),

+                      enrichment_real_norm = (n_real / normFactor_real) / ((n_real / normFactor_real) + n_permuted)) 

+    

+    

+    stopifnot(identical(levels(d2$peak_gene.p.raw.class), levels(d3$peak_gene.p.raw.class)))

+    # 2 enrichment bar plots but combined using facet_wrap

+    d2$set = "r+ / r-"; d3$set = "real / permuted" 

+    d_merged <- tibble::tibble(peak_gene.p.raw.class = c(as.character(d2$peak_gene.p.raw.class), 

+                                                         as.character(d3$peak_gene.p.raw.class)),

+                               ratio = c(d2$ratio_pos_raw, d3$ratio_real_norm),

+                               classAll = c(as.character(d2$class), d3$r_positive),

+                               set = c(d2$set, d3$set)) %>%

+        dplyr::mutate(classAll = factor(classAll, levels=c(paste0("real_",range), paste0("random_",range), "TRUE", "FALSE")),

+                      peak_gene.p.raw.class = factor(peak_gene.p.raw.class, levels = levels(d2$peak_gene.p.raw.class)))

+    

+    d4 = tibble::tibble(peak_gene.p.raw.class = unique(d$peak_gene.p.raw.class), 

+                        n_rpos_real = NA_integer_, n_rpos_random = NA_integer_,

+                        n_rneg_real = NA_integer_, n_rneg_random = NA_integer_,

+                        ratio_permuted_real_rpos_norm = NA_real_,

+                        ratio_permuted_real_rneg_norm = NA_real_)

+    

+    for (i in seq_len(nrow(d4))) {

+        row_d2 = which(d2$class == networkType_details[1] & d2$peak_gene.p.raw.class == d4$peak_gene.p.raw.class[i])

+        stopifnot(length(row_d2) == 1)

+        d4[i, "n_rpos_real"] = d2[row_d2, "sum_pos"] %>% unlist()

+        d4[i, "n_rneg_real"] = d2[row_d2, "sum_neg"] %>% unlist()

+        row_d2 = which(d2$class == paste0("random_",range) & d2$peak_gene.p.raw.class == d4$peak_gene.p.raw.class[i])

+        d4[i, "n_rpos_random"] = d2[row_d2, "sum_pos"] %>% unlist()

+        d4[i, "n_rneg_random"] = d2[row_d2, "sum_neg"] %>% unlist()

+        

+        row_d3 = which(d3$r_positive == TRUE & d3$peak_gene.p.raw.class == d4$peak_gene.p.raw.class[i])

+        d4[i, "ratio_permuted_real_rpos_norm"] = 1- d3[row_d3, "ratio_real_norm"] %>% unlist()

+        row_d3 = which(d3$r_positive == FALSE & d3$peak_gene.p.raw.class == d4$peak_gene.p.raw.class[i])

+        d4[i, "ratio_permuted_real_rneg_norm"] = 1- d3[row_d3, "ratio_real_norm"] %>% unlist()

+    }

+    

+    d4 = d4 %>%

+        dplyr::mutate(ratio_rneg_rpos_real = n_rneg_real / (n_rneg_real + n_rpos_real),

+                      ratio_rneg_rpos_random = n_rneg_random / (n_rneg_random + n_rpos_random),

+                      peak_gene.p.raw.class.bin = as.numeric(peak_gene.p.raw.class)) %>%

+        dplyr::arrange(peak_gene.p.raw.class.bin)

+    

+    d4_melt = reshape2::melt(d4, id  = c("peak_gene.p.raw.class.bin", "peak_gene.p.raw.class")) %>%

+        dplyr::filter(grepl("ratio", variable))

+    

+    

+    list(d = d, d2 = d2, d3 = d3, d4 = d4, d4_melt = d4_melt, d_merged = d_merged)

+    

+}

+

+.classFreq_label <- function(tbl_freq) {

+    paste0(names(tbl_freq), " (", tbl_freq, ")")

+}

+

@@ -58,7 +58,7 @@

     if (!checkmate::testNull(pdfFile)) 

       grDevices::dev.off()

-    futile.logger::flog.info(paste0("Finished writing plots to file ", pdfFile))

+    futile.logger::flog.info(paste0(" Finished writing plots"))

 }

 .startLogger <- function(logfile, level, removeOldLog = TRUE, appenderName = "consoleAndFile", verbose = TRUE) {

@@ -719,3 +719,23 @@ is.installed <- function(mypkg){

     .checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)

   }

 }

+

+

+.getBiomartParameters <- function(genomeAssembly) {

+    

+    .checkOntologyPackageInstallation("biomaRt")

+    

+    if (grepl(x = genomeAssembly, pattern = "^hg\\d\\d")){

+        dataset = "hsapiens_gene_ensembl"

+        if (genomeAssembly == "hg38") {

+            host = "https://www.ensembl.org"

+        } else if (genomeAssembly == "hg19") {

+            host ="https://grch37.ensembl.org"

+        }

+    } else if (grep(x = genomeAssembly, pattern = "^mm\\d\\d")){

+        dataset = "mmusculus_gene_ensembl"

+        host = "https://www.ensembl.org"

+    }

+    

+    list(dataset = dataset, host = host)

+}

\ No newline at end of file

@@ -483,17 +483,13 @@ calculateGeneralEnrichment <- function(GRN, ontology = c("GO_BP", "GO_MF"),

     .checkOntologyPackageInstallation("biomaRt")

-    if (grep(x = GRN@config$parameters$genomeAssembly, pattern = "^hg\\d\\d" )){

-      dataset = "hsapiens_gene_ensembl"

-    } else if (grep(x = GRN@config$parameters$genomeAssembly, pattern = "^mm\\d\\d")){

-      dataset = "mmusculus_gene_ensembl"

-    }

+    params.l = .getBiomartParameters(GRN@config$parameters$genomeAssembly)

     errorOcured = FALSE

     ensembl = tryCatch({ 

-      biomaRt::useEnsembl(biomart = "genes", dataset = dataset)

-      

+        biomaRt::useEnsembl(biomart = "genes", host = params.l[["host"]], dataset = params.l[["dataset"]])

+        

     }, error = function(e) {

       errorOcured = TRUE 

     }

@@ -1019,7 +1015,7 @@ getTopNodes <- function(GRN, nodeType, rankType, n = 0.1, use_TF_gene_network =

 #' 

 #' @inheritParams calculateGeneralEnrichment

 #' @param rankType Character. Default \code{"degree"}. One of: \code{"degree"}, \code{"EV"}, \code{"custom"}. This parameter will determine the criterion to be used to identify the "top" TFs. If set to "degree", the function will select top TFs based on the number of connections to genes they have, i.e. based on their degree-centrality. If set to \code{"EV"} it will select the top TFs based on their eigenvector-centrality score in the network. If set to custom, a set of TF names will have to be passed to the "TF.names" parameter.

-#' @param n Numeric. Default 0.1. If this parameter is passed as a value between (0,1), it is treated as a percentage of top nodes. If the value is passed as an integer it will be treated as the number of top nodes. This parameter is not relevant if \code{rankType = "custom"}.

+#' @param n Numeric. Default 3. If this parameter is passed as a value between 0 and 1, it is treated as a percentage of top nodes. If the value is passed as an integer it will be treated as the number of top nodes. This parameter is not relevant if \code{rankType = "custom"}.

 #' @param TF.names Character vector. Default \code{NULL}. If the rank type is set to \code{"custom"}, a vector of TF names for which the GO enrichment should be calculated should be passed to this parameter.

 #' @return The same \code{\linkS4class{GRN}} object, with the enrichment results stored in the \code{stats$Enrichment$byTF} slot.

 #' @examples 

@@ -1027,7 +1023,7 @@ getTopNodes <- function(GRN, nodeType, rankType, n = 0.1, use_TF_gene_network =

 #' GRN =  loadExampleObject()

 #' GRN =  calculateTFEnrichment(GRN, n = 5, ontology = "GO_BP", forceRerun = FALSE)

 #' @export

-calculateTFEnrichment <- function(GRN, rankType = "degree", n = 0.1, TF.names = NULL,

+calculateTFEnrichment <- function(GRN, rankType = "degree", n = 3, TF.names = NULL,

                                   ontology = c("GO_BP", "GO_MF"), algorithm = "weight01", 

                                   statistic = "fisher", background = "neighborhood",

                                   pAdjustMethod = "BH",

@@ -6,13 +6,15 @@

 #'

 #' @template GRN 

 #' @template outputFolder

-#' @param type Character. Either "peaks" or "rna" or "all". Type of data to plot a PCA for. "peaks" corresponds to the the open chromatin data, while "rna" refers to the RNA-seq counts. If set to "all", PCA will be done for both data modalities. In any case, PCA will be based on the original provided data before any additional normalization has been run (i.e., usually the raw data).

-#' @param topn Integer vector. Default c(500,1000,5000). Number of top variable features to do PCA for. Can be a vector of different numbers (see default).

+#' @param data Character. Either \code{"peaks"} or \code{"rna"} or \code{"all"}. Type of data to plot a PCA for. \code{"peaks"} corresponds to the the open chromatin data, while \code{"rna"} refers to the RNA-seq counts. If set to \code{"all"}, PCA will be done for both data modalities. In any case, PCA will be based on the original provided data before any additional normalization has been run (i.e., usually the raw data).

+#' @param topn Integer vector. Default \code{c(500,1000,5000)}. Number of top variable features to do PCA for. Can be a vector of different numbers (see default).

+#' @param type Character. One of or a combination of \code{"raw"}, \code{"normalized"}, \code{"all"}. Default \code{c("raw", "normalized")}. Should the PCA plots be done based on the raw or normalized data, respectively?

 #' @template forceRerun

 #' @template basenameOutput

 #' @template plotAsPDF

 #' @template pdf_width

 #' @template pdf_height

+#' @template pages

 #' @return The same \code{\linkS4class{GRN}} object, without modifications. In addition, for each specified \code{type}, a PDF file is produced with a PCA. We refer to the Vignettes for details and further explanations.

 #' @examples 

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

@@ -20,17 +22,19 @@

 #' # GRN = plotPCA_all(GRN, topn = 500, type = "rna", plotAsPDF = FALSE)

 #' @export

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

-                        type = c("rna", "peaks"), topn = c(500,1000,5000), 

-                        plotAsPDF = TRUE, pdf_width = 12, pdf_height = 12,

+                        data = c("rna", "peaks"), topn = c(500,1000,5000), type = c("raw", "normalized"), 

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

                         forceRerun = FALSE) {

   GRN = .addFunctionLogToObject(GRN)

   checkmate::assertClass(GRN, "GRN")

-  checkmate::assertSubset(type , c("rna", "peaks", "all"))

+  checkmate::assertSubset(data , c("rna", "peaks", "all"))

+  checkmate::assertSubset(type , c("raw", "normalized", "all"))

   checkmate::assertFlag(plotAsPDF)

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

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

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

   checkmate::assertFlag(forceRerun)

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

@@ -43,7 +47,10 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

   uniqueSamples = GRN@config$sharedSamples

-  if ("rna" %in% type | "all" %in% type) {

+  # Initialize the page counter

+  pageCounter = 1 

+  

+  if ("rna" %in% data | "all" %in% data) {

     # Check if we have raw integer counts. Ust vst transform if, otherwise "regular log2

     transformation = dplyr::if_else(checkmate::testClass(GRN@data$RNA$counts_orig, "DESeqDataSet"), "vst", "log2")

@@ -55,6 +62,10 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

     }

     for (norm in normVector) {

+        

+        if (!(norm %in% type | "all" %in% type)) {

+            next

+        }

       # Only do a transformation if this is raw counts, already normalized counts dont need another transformation

       transformationCur = dplyr::if_else(norm == "raw", transformation, "none")

@@ -85,7 +96,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

                                         " RNA data for all shared samples to file ", fileCur , 

                                         "... This may take a few minutes"))

-        .plot_PCA_wrapper(GRN, matrixCur, transformation = transformationCur, logTransformDensity = logTransformDensity, file = fileCur, topn = topn, pdf_width = pdf_width, pdf_height = pdf_height, plotAsPDF = plotAsPDF)

+        .plot_PCA_wrapper(GRN, matrixCur, transformation = transformationCur, logTransformDensity = logTransformDensity, file = fileCur, topn = topn, pdf_width = pdf_width, pdf_height = pdf_height, pages = pages, plotAsPDF = plotAsPDF)

       } else {

         futile.logger::flog.info(paste0("File ", fileCur, " already exists, not overwriting due to forceRerun = FALSE"))

@@ -95,7 +106,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

   } # end if plot RNA PCA

-  if ("peaks" %in% type | "all" %in% type) {

+  if ("peaks" %in% data | "all" %in% data) {

     # Check if we have raw integer counts

     transformation = dplyr::if_else(checkmate::testClass(GRN@data$peaks$counts_orig, "DESeqDataSet"), "vst", "log2")

@@ -131,7 +142,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

                                         " peaks data for all shared samples to file ", fileCur , 

                                         "... This may take a few minutes"))

-        .plot_PCA_wrapper(GRN, matrixCur, transformation = transformationCur, logTransformDensity = logTransformDensity, file = fileCur, topn = topn, pdf_width = pdf_width, pdf_height = pdf_height, plotAsPDF = plotAsPDF)

+        .plot_PCA_wrapper(GRN, matrixCur, transformation = transformationCur, logTransformDensity = logTransformDensity, file = fileCur, topn = topn, pdf_width = pdf_width, pdf_height = pdf_height, pages = pages, plotAsPDF = plotAsPDF)

       } else {

         futile.logger::flog.info(paste0("File ", fileCur, " already exists, not overwriting due to forceRerun = FALSE"))

@@ -149,7 +160,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

 }

-.plot_PCA_wrapper <- function(GRN, counts, transformation = "vst", scale = TRUE, logTransformDensity, file, topn, pdf_width, pdf_height, plotAsPDF) {

+.plot_PCA_wrapper <- function(GRN, counts, transformation = "vst", scale = TRUE, logTransformDensity, file, topn, pdf_width, pdf_height, pages = NULL, plotAsPDF) {

   checkmate::assertChoice(transformation, c("vst", "log2", "none"))

   checkmate::assertIntegerish(topn, lower = 50, any.missing = FALSE,min.len = 1)

@@ -175,7 +186,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

   futile.logger::flog.info(paste0("Prepare PCA. Count transformation: ", transformation))

-  

+  pageCounter = 1 

   if (plotAsPDF) {

     .checkOutputFile(file)

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

@@ -188,7 +199,11 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

   rv <- matrixStats::rowVars(counts.transf)

   # 0. Density plot for ALL features

-  .plotDensity(counts.transf, logTransform = logTransformDensity, legend = dplyr::if_else(ncol(counts.transf) > 20, FALSE, TRUE))

+  if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+      .plotDensity(counts.transf, logTransform = logTransformDensity, legend = dplyr::if_else(ncol(counts.transf) > 20, FALSE, TRUE))

+      

+  }

+  pageCounter = pageCounter + 1

   for (topn in topn) {

@@ -209,95 +224,105 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

                     pos = seq_len(nPCS))

     # 1. Scree plot

-    g = ggplot(screeplot.df, aes(PCs, variation)) + geom_bar(stat = "identity") + 

-      xlab("Principal component (PC)") + ylab("Explained variation / Cumultative sum (in %)") + 

-      geom_point(data = screeplot.df, aes(pos, variation_sum), color = "red") + 

-      geom_line(data = screeplot.df, aes(pos, variation_sum), color = "red") +

-      ggtitle(paste0("Screeplot for top ", topn, " variable features")) + 

-      theme_bw()

-    

-    plot(g)

-    

+    if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+        

+        g = ggplot(screeplot.df, aes(PCs, variation)) + geom_bar(stat = "identity") + 

+            xlab("Principal component (PC)") + ylab("Explained variation / Cumultative sum (in %)") + 

+            geom_point(data = screeplot.df, aes(pos, variation_sum), color = "red") + 

+            geom_line(data = screeplot.df, aes(pos, variation_sum), color = "red") +

+            ggtitle(paste0("Screeplot for top ", topn, " variable features")) + 

+            theme_bw()

+        

+        plot(g)

+    }

+    pageCounter = pageCounter + 1

     # 2. Metadata correlation with PCs

-    metadata_columns = .findSuitableColumnsToPlot(metadata, remove1LevelFactors = TRUE, verbose = FALSE)

-    .plotPCA_QC(pcaResult = pca, dataCols = colnames(counts.transf), metadataTable = metadata, metadataColumns = metadata_columns, 

-                varn = topn, metadataColumns_fill = metadata_columns, 

-                file = NULL, logTransform = FALSE)

-    

+    if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+        metadata_columns = .findSuitableColumnsToPlot(metadata, remove1LevelFactors = TRUE, verbose = FALSE)

+        .plotPCA_QC(pcaResult = pca, dataCols = colnames(counts.transf), metadataTable = metadata, metadataColumns = metadata_columns, 

+                    varn = topn, metadataColumns_fill = metadata_columns, 

+                    file = NULL, logTransform = FALSE)

+    }

+    pageCounter = pageCounter + 1

     # 3. Correlation plots colored by metadata

     metadata_columns = .findSuitableColumnsToPlot(metadata, remove1LevelFactors = FALSE, verbose = FALSE)

     for (varCur in metadata_columns) {

-      plotTitle = paste0("Top ", topn, " variable features, colored by \n", varCur)

-      

-      skipLegend = FALSE

-      vecCur = metadata[, varCur] %>%unlist()

-      if (is.character(vecCur) & length(unique(vecCur)) > 30) {

-        skipLegend = TRUE

-        plotTitle = paste0(plotTitle, " (legend skipped)")

-      }

-      

-      result = tryCatch({

-        

-        # The following code is taken from plotPCA from DESeq2 and adapted here to be more flexible

-        intgroup.df <- as.data.frame(metadata[, varCur, drop=FALSE])

-        

-        # add the varCur factors together to create a new grouping factor

-        group = metadata[[varCur]]

-        

-        # assembly the data for the plot

-        d <- data.frame(PC1=pca$x[,1], PC2=pca$x[,2], group=group, intgroup.df, name=colnames(counts.transf))

-        

-        d = d %>%

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

-          dplyr::mutate_if(is.logical, as.factor)

-        

-        

-        g = g = ggplot(data=d, aes(x=PC1, y=PC2, color=group)) + geom_point(size=3) + 

-          xlab(paste0("PC1: ",round(percentVar[1] * 100, 1),"% variance")) +

-          ylab(paste0("PC2: ",round(percentVar[2] * 100, 1),"% variance")) +

-          ggtitle(plotTitle) + coord_fixed()

-        

-        

-        

-        if (is.factor(d$group)) {

-          g = g + scale_color_viridis_d()

-          

-        } else {

-          

-          is.date <- function(x) inherits(x, 'Date')

-          

-          if (!is.date(d$group)) {

-            g = g + scale_color_viridis_c()

-          } else {

-            g = g + scale_color_viridis_c(trans = "date")

-          }

-          

-          

+        if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+            

+            plotTitle = paste0("Top ", topn, " variable features, colored by \n", varCur)

+            

+            skipLegend = FALSE

+            vecCur = metadata[, varCur] %>%unlist()

+            if (is.character(vecCur) & length(unique(vecCur)) > 30) {

+                skipLegend = TRUE

+                plotTitle = paste0(plotTitle, " (legend skipped)")

+            }

+            

+            result = tryCatch({

+                

+                # The following code is taken from plotPCA from DESeq2 and adapted here to be more flexible

+                intgroup.df <- as.data.frame(metadata[, varCur, drop=FALSE])

+                

+                # add the varCur factors together to create a new grouping factor

+                group = metadata[[varCur]]

+                

+                # assembly the data for the plot

+                d <- data.frame(PC1=pca$x[,1], PC2=pca$x[,2], group=group, intgroup.df, name=colnames(counts.transf))

+                

+                d = d %>%

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

+                    dplyr::mutate_if(is.logical, as.factor)

+                

+                

+                g = g = ggplot(data=d, aes(x=PC1, y=PC2, color=group)) + geom_point(size=3) + 

+                    xlab(paste0("PC1: ",round(percentVar[1] * 100, 1),"% variance")) +

+                    ylab(paste0("PC2: ",round(percentVar[2] * 100, 1),"% variance")) +

+                    ggtitle(plotTitle) + coord_fixed()

+                

+                

+                

+                if (is.factor(d$group)) {

+                    g = g + scale_color_viridis_d()

+                    

+                } else {

+                    

+                    is.date <- function(x) inherits(x, 'Date')

+                    

+                    if (!is.date(d$group)) {

+                        g = g + scale_color_viridis_c()

+                    } else {

+                        g = g + scale_color_viridis_c(trans = "date")

+                    }

+                    

+                    

+                }

+                

+                

+                if (skipLegend) {

+                    g = g + theme(legend.position = "none")

+                } 

+                

+                plot(g)

+                

+            }, error = function(e) {

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

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

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

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

+            })

         }

-        

-        

-        if (skipLegend) {

-          g = g + theme(legend.position = "none")

-        } 

-        

-        plot(g)

-        

-      }, error = function(e) {

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

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

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

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

-      })

+        pageCounter = pageCounter + 1

-    }

+    } # end for each variable

   }

+  .printWarningPageNumber(pages, pageCounter)

   if (plotAsPDF) grDevices::dev.off()

 }

@@ -426,7 +451,7 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

     cell_fun = function(j, i, x, y, width, height, fill) {

       grid::grid.text(sprintf("%.2f", r2_cv[i, j]), x, y, gp = grid::gpar(fontsize = 10))

     }

-  )

+  ) %>% plot()

@@ -450,8 +475,12 @@ plotPCA_all <- function(GRN, outputFolder = NULL, basenameOutput = NULL,

 #' @template outputFolder

 #' @template basenameOutput

 #' @template plotDetails

-#' @param plotPermuted \code{TRUE} or \code{FALSE}. Default  \code{TRUE}. Also produce the diagnostic plots for permuted data?

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

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

+#' @template plotAsPDF

+#' @template pdf_width

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

+#' @template pages

 #' @template forceRerun

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

 #' @examples 

@@ -463,8 +492,9 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

                                         outputFolder = NULL, 

                                         basenameOutput = NULL, 

                                         plotDetails = FALSE,

-                                        plotPermuted = TRUE,

+                                        dataType = c("real", "permuted"),

                                         nTFMax = NULL,

+                                        plotAsPDF = TRUE, pdf_width = 12, pdf_height_base = 8, pages = NULL,

                                         forceRerun = FALSE) {

   GRN = .addFunctionLogToObject(GRN)

@@ -473,8 +503,12 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

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

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

   checkmate::assertFlag(plotDetails)

-  checkmate::assertFlag(plotPermuted)

+  checkmate::assertSubset(dataType, c("real", "permuted"))

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

+  checkmate::assertFlag(plotAsPDF)

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

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

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

   checkmate::assertFlag(forceRerun)

   useGCCorrection = GRN@config$parameters$useGCCorrection

@@ -482,29 +516,38 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

   outputFolder = .checkOutputFolder(GRN, outputFolder)

+  dataType2 = c()

+  if ("real" %in% dataType) dataType2 = c(0)

+  if ("permuted" %in% dataType) dataType2 = c(dataType2, 1)

+  

   for (permutationCur in 0:.getMaxPermutation(GRN)) {

-    if (!plotPermuted & permutationCur != 0) {

+    if (!permutationCur %in% dataType2) {

       next

     }

     suffixFile = .getPermutationSuffixStr(permutationCur)

-    fileCur = paste0(outputFolder, 

-                     dplyr::if_else(is.null(basenameOutput), .getOutputFileName("plot_TFPeak_fdr"), basenameOutput), 

+    fileCur = paste0(outputFolder, dplyr::if_else(is.null(basenameOutput), .getOutputFileName("plot_TFPeak_fdr"), basenameOutput), 

                      suffixFile, ".pdf")

-    if (!file.exists(fileCur) | forceRerun) {

+    

+    if (!file.exists(fileCur) | !plotAsPDF | forceRerun) {

+        

+        if (!plotAsPDF) fileCur = NULL

-      heightCur = 8* length(GRN@config$TF_peak_connectionTypes)

-      .plot_TF_peak_fdr(GRN, perm = permutationCur, useGCCorrection = useGCCorrection, 

-                        plotDetails = plotDetails, fileCur, width = 7, height = heightCur,

-                        nPagesMax = nTFMax) 

+        heightCur = pdf_height_base * length(GRN@config$TF_peak_connectionTypes)

+        .plot_TF_peak_fdr(GRN, perm = permutationCur, useGCCorrection = useGCCorrection, 

+                        plotDetails = plotDetails, fileCur, width = pdf_width, height = heightCur,

+                        nPagesMax = nTFMax, pages = pages) 

     }

+    # TODO: page selection not implemented here yet

     fileCur = paste0(outputFolder, .getOutputFileName("plot_TFPeak_fdr_GC"), suffixFile)

-    if (useGCCorrection & (!file.exists(fileCur) | forceRerun)) {

-      .plotTF_peak_GC_diagnosticPlots(GRN, perm = permutationCur, fileCur, width = 7, height = 8) 

+    if (useGCCorrection & (!file.exists(fileCur) | !plotAsPDF | forceRerun)) {

+        

+        if (!plotAsPDF) fileCur = NULL

+        .plotTF_peak_GC_diagnosticPlots(GRN, perm = permutationCur, fileCur, width = pdf_width, height = pdf_height_base) 

     }

     # TODO: handle multiple activities and actually write this function

@@ -519,7 +562,7 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

 }

-

+# TODO: Currently not used anywhere

 .generateTF_GC_diagnosticPlots <- function(TFCur, GC_classes_foreground.df, GC_classes_background.df, GC_classes_all.df, peaksForeground, peaksBackground, peaksBackgroundGC) {

   GC_classes_background_GC.df = peaksBackgroundGC %>%

@@ -627,15 +670,15 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

 }

-.plot_TF_peak_fdr <- function(GRN, perm, useGCCorrection, plotDetails = FALSE, file = NULL, width = 7, height = 7, nPagesMax = NULL) {

+.plot_TF_peak_fdr <- function(GRN, perm, useGCCorrection, plotDetails = FALSE, file = NULL, width = 7, height = 7, nPagesMax = NULL, pages = NULL) {

   start = Sys.time()

   futile.logger::flog.info(paste0("Plotting FDR curves for each TF", dplyr::if_else(is.null(file), "", paste0(" to file ", file))))

-  if (!is.null(file)) {

-    plots.l = list()

-    index = 1

-  } 

+  pageCounter = 1  

+  

+  plots.l = list()

+  index = 1

   # Dont take all TF, some might be missing.

   connections_TF_peak = GRN@connections$TF_peaks[[as.character(perm)]]$connectionStats

@@ -755,35 +798,42 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

       }

-      

-      

-      if (is.null(file)) {

-        plot(g)

-      } else {

-        plots.l[[index]] = g

-        index = index + 1

-        

-        if (plotDetails) {

-          plots.l[[index]] = g2

+      # We now created either just one plot g or g + g2 for "one" direction

+      if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+          plots.l[[index]] = g

           index = index + 1

-        }

-        

       }

+

+      if (plotDetails) {

+          

+          if (is.null(pages) | (!is.null(pages) && pageCounter %in% pages)) {

+              plots.l[[index]] = g2

+              index = index + 1

+          }

+

+          # Only here we already have a new page, as we produced 2 plots in this iteration alone

+          pageCounter = pageCounter + 1 

+       }

+

+    } # end for both directions

+    

+    if (!plotDetails) {

+        pageCounter = pageCounter + 1 

     }

+

+    

   } # end allTF

-  if (!is.null(file)) {

-    futile.logger::flog.info(paste0(" Finished generating plots, start plotting to file ",file, ". This may take a few minutes."))

-    .printMultipleGraphsPerPage(plots.l, nCol = 1, nRow = 2, pdfFile = file, height = height, width = width)

-    

-  }

+  .printWarningPageNumber(pages, pageCounter)

+  

+  .printMultipleGraphsPerPage(plots.l, nCol = 1, nRow = 2, pdfFile = file, height = height, width = width)

   .printExecutionTime(start)

 }

-.plotTF_peak_TFActivity_QC <- function(GRN, perm = permIndex, fileCur, width = 7, height = 8) {

+.plotTF_peak_TFActivity_QC <- function(GRN, perm = 0, fileCur, width = 7, height = 8) {

   # TODO

 }

@@ -796,13 +846,14 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

 #' @template GRN 

 #' @template outputFolder

 #' @template basenameOutput

-#' @param gene.types List of character vectors. Default list(c("all"), c("protein_coding"), c("protein_coding", "lincRNA")). Vectors of gene types to consider for the diagnostic plots. Multiple distinct combinations of gene types can be specified. For example, if set to \code{list(c("protein_coding", "lincRNA"), c("protein_coding"), c("all"))}, 3 distinct PDFs will be produced, one for each element of the list. The first file would only consider protein-coding and lincRNA genes, while the second plot only considers protein-coding ones. The special keyword "all" denotes all gene types found (usually, there are many gene types present, also more exotic and rare ones).

+#' @param gene.types List of character vectors. Default list(c("protein_coding", "lincRNA")). Vectors of gene types to consider for the diagnostic plots. Multiple distinct combinations of gene types can be specified. For example, if set to \code{list(c("protein_coding", "lincRNA"), c("protein_coding"), c("all"))}, 3 distinct PDFs will be produced, one for each element of the list. The first file would only consider protein-coding and lincRNA genes, while the second plot only considers protein-coding ones. The special keyword "all" denotes all gene types found (usually, there are many gene types present, also more exotic and rare ones).

 #' @param useFiltered Logical. TRUE or FALSE. Default FALSE. If set to \code{FALSE}, the diagnostic plots will be produced based on all peak-gene connections. This is the default and will usually be best to judge whether the background behaves as expected. If set to TRUE, the diagnostic plots will be produced based on the filtered set of connections. For this, the function \code{link{filterGRNAndConnectGenes}} must have been run before.

 #' @template plotDetails

 #' @param plotPerTF Logical. TRUE or FALSE. Default \code{FALSE}. If set to \code{FALSE}, the diagnostic plots will be done across all TF (the default), while setting it to \code{TRUE} will generate the QC plots TF-specifically, including "all" TF, sorted by the number of connections.

 #' @template plotAsPDF

 #' @template pdf_width

 #' @template pdf_height

+#' @template pages

 #' @template forceRerun

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

 #' @examples 

@@ -815,11 +866,11 @@ plotDiagnosticPlots_TFPeaks <- function(GRN,

 plotDiagnosticPlots_peakGene <- function(GRN, 

                                          outputFolder = NULL, 

                                          basenameOutput = NULL, 

-                                         gene.types = list(c("all"), c("protein_coding"), c("protein_coding", "lincRNA")), 

+                                         gene.types = list(c("protein_coding", "lincRNA")), 

                                          useFiltered = FALSE, 

                                          plotDetails = FALSE,

                                          plotPerTF = FALSE,

-                                         plotAsPDF = TRUE, pdf_width = 12, pdf_height = 12,

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

                                          forceRerun= FALSE) {

   GRN = .addFunctionLogToObject(GRN)

@@ -834,6 +885,7 @@ plotDiagnosticPlots_peakGene <- function(GRN,

   checkmate::assertFlag(plotAsPDF)

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

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

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

   checkmate::assertFlag(forceRerun)

@@ -864,20 +916,9 @@ plotDiagnosticPlots_peakGene <- function(GRN,

                                     plotPerTF = plotPerTF,

                                     pdf_width = pdf_width, pdf_height = pdf_height,

                                     plotDetails = plotDetails,

+                                    pages = pages,

                                     forceRerun = forceRerun)

-  # Summarize all data, both real and random

-  # filteredStr = dplyr::if_else(useFiltered, "_filtered", "")

-  # for (geneTypesSelected in gene.types) {

-  #   

-  #   filenameCur = paste0(outputFolder, .getOutputFileName("plot_peakGene_diag"), "_", paste0(geneTypesSelected, collapse = "+"), filteredStr, ".pdf")

-  #   if (!file.exists(filenameCur) | forceRerun) {

-  #    

-  #   } else {

-  #     futile.logger::flog.info(paste0("File ", filenameCur, " already exists. As forceRerun = FALSE, not overwriting it."))

-  #   }

-  #   

-  # }

   GRN

@@ -891,6 +932,7 @@ plotDiagnosticPlots_peakGene <- function(GRN,

                                               fileBase = NULL,

                                               pdf_width = 12,

                                               pdf_height = 12,

+                                              pages = NULL,

                                               forceRerun = FALSE) {

   start = Sys.time()

@@ -1005,8 +1047,7 @@ plotDiagnosticPlots_peakGene <- function(GRN,

     freq_class = gsub(freq_class, pattern = "random", replacement = "permuted")

     names(freq_class) <- names(freqs)

-    

-    

+

     xlabels_peakGene_r.class = levels(peakGeneCorrelations.all$peak_gene.r.class)

     nCur = length(xlabels_peakGene_r.class)

     xlabels_peakGene_r.class[setdiff(seq_len(nCur), c(1, floor(nCur/2), nCur))] <- ""

@@ -1020,9 +1061,14 @@ plotDiagnosticPlots_peakGene <- function(GRN,

     nCur = length(xlabels_peakGene_praw.class)

     xlabels_peakGene_praw.class[setdiff(seq_len(nCur), c(1, floor(nCur/2), nCur))] <- ""

-    

+    #

+    # ITERATE THROUGH ALL GENE TYPES, WITH ONE PER PLOT

+    #

     for (geneTypesSelected in gene.types) {

+      # Reset page counter for each PDF anew

+      pageCounter = 1

+        

       futile.logger::flog.info(paste0(" Gene type ", paste0(geneTypesSelected, collapse = "+")))

       if ("all" %in% geneTypesSelected) {

@@ -1064,116 +1110,120 @@ plotDiagnosticPlots_peakGene <- function(GRN,