Commit 46e702ac authored by Christian Arnold's avatar Christian Arnold
Browse files

GRaNIE: New version 0.14

parent 214f1b28
Pipeline #28081 passed with stage
in 13 seconds
^.*\.Rproj$
^\.Rproj\.user$
^doc$
^Meta$
^_pkgdown\.yml$
^docs$
^pkgdown$
Package: GRN
Title: GRN: Reconstruction cell type specific gene regulatory networks using chromatin accessibility and RNA-seq data
Version: 0.10.1
Package: GRaNIE
Title: GRaNIE: Reconstruction cell type specific gene regulatory networks including enhancers using chromatin accessibility and RNA-seq data
Version: 0.14.0
Encoding: UTF-8
Authors@R: c(person("Christian", "Arnold", email =
"christian.arnold@embl.de", role = c("cre","aut")),
person("Judith", "Zaugg", email =
"judith.zaugg@embl.de", role = "aut"))
Author: Christian Arnold [aut, cre], Judith Zaugg [aut]
Maintainer: Christian Arnold <christian.arnold@embl.de>
"judith.zaugg@embl.de", role = c("aut")),
person("Rim", "Moussa", email =
"rim.moussa01@gmail.com", role = "ctb"),
person("Armando", "Reyes-Palomares", email =
"armandorp@gmail.com", role = "ctb"),
person("Giovanni", "Palla", email =
"giov.pll@gmail.com", role = "ctb"),
person("Maksim", "Kholmatov", email =
"maksim.kholmatov@embl.de", role = "ctb"))
Description: Genetic variants associated with diseases often affect non-coding regions, thus likely having a regulatory role. To understand the effects of genetic variants in these regulatory regions, identifying genes that are modulated by specific regulatory elements (REs) is crucial. The effect of gene regulatory elements, such as enhancers, is often cell-type specific, likely because the combinations of transcription factors (TFs) that are regulating a given enhancer have celltype specific activity. This TF activity can be quantified with existing tools such as diffTF and captures differences in binding of a TF in open chromatin regions. Collectively, this forms a gene regulatory network (GRN) with cell-type and data-specific TF-RE and RE-gene links. Here, we reconstruct such a GRN using bulk RNAseq and open chromatin (e.g., using ATACseq or ChIPseq for open chromatin marks) and optionally TF activity data. Our network contains different types of links, connecting TFs to regulatory elements, the latter of which is connected to genes in the vicinity or within the same chromatin domain (TAD). We use a statistical framework to assign empirical FDRs and weights to all links using a permutation-based approach.
Imports:
futile.logger,
......@@ -21,6 +27,7 @@ Imports:
RColorBrewer,
colorspace,
pheatmap,
ComplexHeatmap,
limma,
DESeq2,
csaw,
......@@ -29,7 +36,6 @@ Imports:
utils,
methods,
stringr,
dplyr,
scales,
BiocManager,
BiocParallel,
......@@ -42,22 +48,24 @@ Imports:
IRanges,
SummarizedExperiment,
forcats,
ggplotify,
gridExtra,
limma,
purrr,
tibble,
tidyr,
tidyselect,
readr,
grid
Depends:
R (>= 3.6),
tidyverse,
topGO,
tidyr,
dplyr,
stats,
grDevices,
graphics,
geneplotter,
readr,
vsn
Depends:
R (>= 3.5),
tidyverse
magrittr,
tibble
Suggests:
knitr,
BSgenome.Hsapiens.UCSC.hg19,
......@@ -68,11 +76,24 @@ Suggests:
TxDb.Hsapiens.UCSC.hg38.knownGene,
TxDb.Mmusculus.UCSC.mm10.knownGene,
TxDb.Mmusculus.UCSC.mm9.knownGene,
IHW
IHW,
ChIPseeker,
testthat (>= 3.0.0)
VignetteBuilder: knitr
biocViews: Software
biocViews:
Software,
GeneExpression,
GeneRegulation,
NetworkInference,
GeneSetEnrichment,
BiomedicalInformatics,
Genetics,
Transcriptomics,
ATACSeq,
RNASeq
License: LGPL (>= 3)
LazyData: true
URL: https://grp-zaugg.embl-community.io/grn
BugReports: https://git.embl.de/grp-zaugg/grn/issues
URL: https://grp-zaugg.embl-community.io/GRaNIE
BugReports: https://git.embl.de/grp-zaugg/GRaNIE/issues
RoxygenNote: 7.1.1
Config/testthat/edition: 3
......@@ -4,8 +4,13 @@ export(AR_classification_wrapper)
export(addConnections_TF_peak)
export(addConnections_peak_gene)
export(addData)
export(addData_TFActivity)
export(addTFBS)
export(add_TF_gene_correlation)
export(calculateCommunitiesEnrichment)
export(calculateCommunitiesStats)
export(calculateGeneralEnrichment)
export(calculateTFEnrichment)
export(deleteIntermediateData)
export(filterData)
export(filterGRNAndConnectGenes)
......@@ -13,22 +18,33 @@ export(generateStatsSummary)
export(getCounts)
export(getGRNConnections)
export(getParameters)
export(getTopNodes)
export(importTFData)
export(initializeGRN)
export(nGenes)
export(nPeaks)
export(overlapPeaksAndTFBS)
export(performAllNetworkAnalyses)
export(plotCommunitiesEnrichment)
export(plotCommunitiesStats)
export(plotDiagnosticPlots_TFPeaks)
export(plotDiagnosticPlots_peakGene)
export(plotGeneralEnrichment)
export(plotGeneralGraphStats)
export(plotPCA_all)
export(plotTFEnrichment)
export(plot_stats_connectionSummary)
import(GenomicRanges)
import(IRanges)
import(checkmate)
import(ggplot2)
import(graphics)
import(grid)
import(patchwork)
importFrom(BiocParallel,MulticoreParam)
importFrom(BiocParallel,bplapply)
importFrom(BiocParallel,bpok)
importFrom(BiocParallel,multicoreWorkers)
import(rlang)
import(tibble)
import(topGO)
import(utils)
importFrom(GenomicRanges,mcols)
importFrom(methods,new)
importFrom(rlang,.data)
......
<!--
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-->
## GRN 0.10 (2021-06-01)
## GRaNIE 0.9 - 0.14 (2021-12-13)
### Major changes
- decreased memory footprint and increased speed for *plotDiagnosticPlots_peakGene*
- major overhaul and continuous work on peak-gene QC plots
- the *filterData* functions has now more filter parameters, such as filtering for CV. Also, all filters uniformly have a *min* and *max* filter.
- integrated network statistics and various enrichment analyses
- handling of edge cases and rare events in various functions
- packages have been renamed to *GRaNIE* as basename (before: *GRN*)
### Bug fixes
- various minor bug fixes (e.g. genes with identical identifiers are now correctly summarized)
- various minor bug fixes
### Minor changes
- replaced the function argument `permutation` to `permuted` ( TRUE or FALSE) in all exported functions, which is easier to understand for users
- changed the object structure slightly and moved some gene and peak annotation data (such as mean, CV) to the appropriate annotation slot
## GRN 0.8 (2021-05-07)
## GRaNIE 0.8 (2021-05-07)
### Major changes
- improved PCA plotting, PCA plots are now produced for both raw and normalized data
- new filters for the function *filterGRNAndConnectGenes* (*peak_gene.maxDistance*) as well as more flexibility how to adjust the peak-gene raw p-values for multiple testing (including the possibility to use IHW - experimental)
- new filters for the function *filterGRaNIEAndConnectGenes* (*peak_gene.maxDistance*) as well as more flexibility how to adjust the peak-gene raw p-values for multiple testing (including the possibility to use IHW - experimental)
- new function *plotDiagnosticPlots_TFPeaks* for plotting (this function was previously called only internally, but is now properly exported), in analogy to *plotDiagnosticPlots_peakGene*
### Bug fixes
......@@ -36,7 +35,7 @@
- some functions have been added / renamed to make the workflow more clear and streamlined, see Vignette for details
- some default parameters changed
## GRN 0.7 (2021-03-12)
## GRaNIE 0.7 (2021-03-12)
### Major changes
......@@ -51,10 +50,8 @@
### Minor changes
<!--
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-->
## GRN 0.6 (2021-02-09)
## GRaNIE 0.6 (2021-02-09)
### Major changes
......@@ -67,9 +64,7 @@
### Minor changes
<!--
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-->
## GRN 0.5 (2021-02-02)
## GRaNIE 0.5 (2021-02-02)
first published package version
#### S4 class definition
#' Create, represent, investigate, quantify and visualize gene regulatory networks (GRNs)
#' Create, represent, investigate, quantify and visualize enhancer-mediated gene regulatory networks (\strong{eGRNs})
#'
#' The class \code{\linkS4class{GRN}} stores data and information related to our GRN approach to construct gene regulatory networks out of open chromatin and RNA-Seq data. See the description below for more details, and visit our project website at https://grp-zaugg.embl-community.io/grn and have a look at the various Vignettes.
#' The class \code{\linkS4class{GRN}} stores data and information related to our \code{eGRN} approach to construct enhancer-mediated gene regulatory networks out of open chromatin and RNA-Seq data. See the description below for more details, and \strong{visit our project website at https://grp-zaugg.embl-community.io/GRaNIE and have a look at the various Vignettes}.
#' @section Constructors:
#' Currently, a \code{\linkS4class{GRN}} object is created by executing the function \code{\link{initializeGRN}}.
#'
#' @section Accessors:
#' In the following code snippets, \code{GRN} is a \code{\linkS4class{GRN}} object. \strong{In addition, see the workflow vignette (browseVignettes("GRN") for a full workflow that uses all accessors.}
#' In the following code snippets, \code{GRN} is a \code{\linkS4class{GRN}} object.
#'
#'# Get general annotation of a GRN object
#'
......@@ -101,7 +101,7 @@ setValidity("GRN", .validGRNObj)
# isDev = TRUE
# )
new(getClassDef("GRN", package = utils::packageName(), inherits = FALSE))
new(methods::getClassDef("GRN", package = utils::packageName(), inherits = FALSE))
}
# #
......@@ -1145,7 +1145,7 @@ setMethod("show",
#' @return Integer. Number of peaks hat are defined in the \code{\linkS4class{GRN}} object, either by excluding (filter = TRUE) or including (filter = FALSE) peaks that are currently marked as \emph{filtered} (see method TODO)
#' @examples
# #' data(GRN, package="GRN")
# #' nPeaks(GRN, filter = TRUE)
#' nPeaks(GRN, filter = TRUE)
# #' nPeaks(GRN, filter = FALSE)
#' @export
#' @aliases peaks
......@@ -1179,7 +1179,7 @@ nPeaks <- function(GRN, filter = TRUE) {
#' @return Integer. Number of genes hat are defined in the \code{\linkS4class{GRN}} object, either by excluding (filter = TRUE) or including (filter = FALSE) genes that are currently marked as \emph{filtered} (see method TODO)
#' @examples
# #' data(GRN, package="GRN")
# #' nGenes(GRN, filter = TRUE)
#' nGenes(GRN, filter = TRUE)
# #' nGenes(GRN, filter = FALSE)
#' @export
#' @aliases genes
......
#' GRN: Reconstruction cell type specific gene regulatory networks using chromatin accessibility and RNA-seq data (general package information)
#' \strong{GRaNIE} (\strong{G}ene \strong{R}egul\strong{a}tory \strong{N}etwork \strong{I}nference including \strong{E}nhancers): Reconstruction and evaluation of data-driven, cell type specific gene regulatory networks including enhancers using chromatin accessibility and RNAseq data (general package information)
#'
#' Genetic variants associated with diseases often affect non-coding regions, thus likely having a regulatory role. To understand the effects of genetic variants in these regulatory regions, identifying genes that are modulated by specific regulatory elements (REs) is crucial. The effect of gene regulatory elements, such as enhancers, is often cell-type specific, likely because the combinations of transcription factors (TFs) that are regulating a given enhancer have celltype specific activity. This TF activity can be quantified with existing tools such as diffTF and captures differences in binding of a TF in open chromatin regions. Collectively, this forms a gene regulatory network (GRN) with cell-type and data-specific TF-RE and RE-gene links. Here, we reconstruct such a GRN using bulk RNAseq and open chromatin (e.g., using ATACseq or ChIPseq for open chromatin marks) and optionally TF activity data. Our network contains different types of links, connecting TFs to regulatory elements, the latter of which is connected to genes in the vicinity or within the same chromatin domain (TAD). We use a statistical framework to assign empirical FDRs and weights to all links using a permutation-based approach.
#' Genetic variants associated with diseases often affect non-coding regions, thus likely having a regulatory role. To understand the effects of genetic variants in these regulatory regions, identifying genes that are modulated by specific regulatory elements (REs) is crucial. The effect of gene regulatory elements, such as enhancers, is often cell-type specific, likely because the combinations of transcription factors (TFs) that are regulating a given enhancer have celltype specific activity. This TF activity can be quantified with existing tools such as \code{diffTF} and captures differences in binding of a TF in open chromatin regions. Collectively, this forms a gene regulatory network (\code{eGRN}) with cell-type and data-specific TF-RE and RE-gene links. Here, we reconstruct such a \code{eGRN} using bulk RNAseq and open chromatin (e.g., using ATACseq or ChIPseq for open chromatin marks) and optionally TF activity data. Our network contains different types of links, connecting TFs to regulatory elements, the latter of which is connected to genes in the vicinity or within the same chromatin domain (TAD). We use a statistical framework to assign empirical FDRs and weights to all links using a permutation-based approach.
#'
#' @section Package functions:
#' See the Vignettes for a workflow example and more generally \url{https://grp-zaugg.embl-community.io/grn/articles/} for all project-related information.
#' See the Vignettes for a workflow example and more generally \url{https://grp-zaugg.embl-community.io/GRaNIE/articles/} for all project-related information.
#'
#' @section GRN object:
#' The GRN* packages work with GRN objects. See \code{\linkS4class{GRN}} for details.
#' The \code{GRaNIE} package works with \code{GRN} objects. See \code{\linkS4class{GRN}} for details.
#'
#' @section Contact Information:
#' Please check out \url{https://grp-zaugg.embl-community.io/grn} for how to get in contact with us.
#' Please check out \url{https://grp-zaugg.embl-community.io/GRaNIE} for how to get in contact with us.
#'
#' @docType package
#' @keywords GRN, GRN-package
#' @name GRN
#' @keywords GRaNIE, GRaNIE-package
#' @name GRaNIE
NULL
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -101,7 +101,7 @@
dd <- DESeq2::DESeqDataSetFromMatrix(countData = countDataNew,
colData = sampleData.df,
design = as.formula(" ~ 1"))
design = stats::as.formula(" ~ 1"))
dd = DESeq2::estimateSizeFactors(dd)
counts.norm = DESeq2::counts(dd, normalized = TRUE)
......@@ -248,6 +248,11 @@
start = Sys.time()
futile.logger::flog.info(paste0("Calculate classification thresholds for repressors / activators"))
act.rep.thres.l = list()
if (is.null(par.l$internal$allClassificationThresholds)) {
message = paste0("GRN object needs to be updated. Slot internal$allClassificationThresholds is NULL.")
.checkAndLogWarningsAndErrors(NULL, message, isWarning = FALSE)
}
for (thresholdCur in par.l$internal$allClassificationThresholds) {
act.rep.cur = quantile(sort(apply(background, MARGIN = 2, FUN = .my.median)),
......@@ -330,8 +335,8 @@
}
# Removing NAs actually makes a difference, as these are "artifical" anyway here due to the two matrices let's remove them
dataMotif = na.omit(t.cor.sel.matrix[,TFCur])
dataBackground = na.omit(t.cor.sel.matrix.non[,TFCur])
dataMotif = stats::na.omit(t.cor.sel.matrix[,TFCur])
dataBackground = stats::na.omit(t.cor.sel.matrix.non[,TFCur])
# Skip if NA for median correlation
if (is.na(output.global.TFs[rowNo, colnameMedianCor])) {
......@@ -346,7 +351,7 @@
alternativeTest = "less"
}
testResults = wilcox.test(dataMotif, dataBackground, alternative = alternativeTest)
testResults = stats::wilcox.test(dataMotif, dataBackground, alternative = alternativeTest)
stopifnot(length(rowNo) == 1)
output.global.TFs[rowNo,colnameClassificationPVal] = testResults$p.value
......@@ -393,7 +398,7 @@
colnamesIndex = which(grepl("final", colnames(output.global.TFs)))
for (colnameCur in colnamesIndex) {
futile.logger::flog.info(paste0(" Column ", colnames(output.global.TFs)[colnameCur]))
tbl = table(output.global.TFs %>% pull(colnameCur))
tbl = table(output.global.TFs %>% dplyr::pull(colnameCur))
futile.logger::flog.info(paste0(" ", paste0(names(tbl), ": ", tbl), collapse = ", "))
}
......@@ -408,7 +413,7 @@
##################
# Density plots for TFs
#' @import graphics
.plot_density <- function(foreground.m, background.m, corMethod, file = NULL, ...) {
start = Sys.time()
......@@ -421,7 +426,7 @@
n_notNA = length(which(!is.na(foreground.m[,colCur])))
if (n_notNA > 1){
yMaxCur = max(c(density(foreground.m[,colCur], na.rm = TRUE)$y, density(background.m[,colCur], na.rm = T)$y))
yMaxCur = max(c(stats::density(foreground.m[,colCur], na.rm = TRUE)$y, stats::density(background.m[,colCur], na.rm = TRUE)$y))
if (yMaxCur > yMax) {
yMax = yMaxCur + 0.1
......@@ -431,7 +436,7 @@
}
if (!is.null(file)) {
pdf(file, ...)
grDevices::pdf(file, ...)
}
for (colCur in seq_len(ncol(foreground.m))) {
......@@ -449,14 +454,14 @@
# Show n for both foreground and background
legendAnno = c(paste0("Motif (n=", .prettyNum(nForeground), ")"), paste0("Non-motif (n=", .prettyNum(nBackground), ")"))
plot(density(dataMotif, na.rm = TRUE), xlim = c(-1,1), ylim = c(0,yMax),
plot(stats::density(dataMotif, na.rm = TRUE), xlim = c(-1,1), ylim = c(0,yMax),
main= mainLabel, lwd=2.5, col="red", axes = FALSE, xlab = paste0(.firstLetterUppercase(corMethod), " correlation"))
abline(v=0, col="black", lty=2)
legend("topleft",box.col = adjustcolor("white",alpha.f = 0), legend = legendAnno, lwd = c(2,2),cex = 0.8, col = c("red","darkgrey"), lty = c(1,1) )
legend("topleft",box.col = grDevices::adjustcolor("white",alpha.f = 0), legend = legendAnno, lwd = c(2,2),cex = 0.8, col = c("red","darkgrey"), lty = c(1,1) )
axis(side = 1, lwd = 1, line = 0)
axis(side = 2, lwd = 1, line = 0, las = 1)
lines(density(dataBackground, na.rm = T), lwd = 2.5, col = "darkgrey")
lines(stats::density(dataBackground, na.rm = TRUE), lwd = 2.5, col = "darkgrey")
} else {
futile.logger::flog.warn(paste0(" Not enough data for estimating densities for TF ", TFCur, ", skip for plotting."))
......@@ -466,7 +471,7 @@
}
if (!is.null(file)) {
dev.off()
grDevices::dev.off()
}
.printExecutionTime(start)
......@@ -481,7 +486,7 @@
if (!is.null(file)) {
pdf(file, ...)
grDevices::pdf(file, ...)
}
DESeq_results = DESeq2::results(DESeq_obj) %>%
......@@ -515,7 +520,7 @@
cor.res.l = list()
for (corMethodCur in c("pearson", "spearman")) {
cor.res.l[[corMethodCur]] = cor.test(output.global.TFs.cur$weighted_meanDifference,
cor.res.l[[corMethodCur]] = stats::cor.test(output.global.TFs.cur$weighted_meanDifference,
output.global.TFs.cur$log2FoldChange,
method = corMethodCur)
}
......@@ -537,14 +542,14 @@
}
if (!is.null(file)) {
dev.off()
grDevices::dev.off()
}
.printExecutionTime(start)
}
#' @import graphics
.plot_AR_thresholds <- function(median.cor.tfs, median.cor.tfs.non, par.l, act.rep.thres.l, corMethod, file = NULL, ...) {
start = Sys.time()
......@@ -563,7 +568,7 @@
if (!is.null(file)) {
pdf(file, ...)
grDevices::pdf(file, ...)
}
par(mfrow = c(1,1))
......@@ -576,10 +581,10 @@
mainCur = paste0("Stringency: ", thresCur)
plot(median.cor.tfs.non, 1:nTF,
plot(median.cor.tfs.non, seq_len(nTF),
xlim = xlim, ylim = ylim, main=mainCur, xlab=xlab, ylab=ylab,
col = adjustcolor("darkgrey",alpha=1), pch = 16, cex = 0.5,axes = FALSE)
points(median.cor.tfs, 1:nTF,
col = grDevices::adjustcolor("darkgrey",alpha=1), pch = 16, cex = 0.5,axes = FALSE)
points(median.cor.tfs, seq_len(nTF),
pch=16, cex = 0.5,
col= dplyr::case_when(median.cor.tfs > thresCur.v[2] ~ par.l$internal$colorCategories["activator"],
median.cor.tfs < thresCur.v[1] ~ par.l$internal$colorCategories["repressor"],
......@@ -617,7 +622,7 @@
}
if (!is.null(file)) {
dev.off()
grDevices::dev.off()
}
.printExecutionTime(start)
......@@ -639,7 +644,7 @@
}
if (!is.null(file)) {
pdf(file, ...)
grDevices::pdf(file, ...)
}
cor.r.filt.m <- sort.cor.m[,as.character(HOCOMOCO_mapping.df.exp$HOCOID)]
......@@ -653,7 +658,7 @@
TF_Peak_all.m <- TF.peakMatrix.df
TF_Peak.m <- TF_Peak_all.m
for (i in 1:ncol( cor.r.filt.m)) {
for (i in seq_len(ncol(cor.r.filt.m))) {
TF = colnames( cor.r.filt.m)[i]
## for the background, use all peaks
h_noMotif = hist( cor.r.filt.m[,TF][TF_Peak_all.m[,TF] == 0], breaks = BREAKS, plot = FALSE)
......@@ -744,7 +749,7 @@
if (!is.null(file)) {
dev.off()
grDevices::dev.off()
}
.printExecutionTime(start)
......
......@@ -48,14 +48,14 @@
if (!checkmate::testNull(pdfFile)) {
.clearOpenDevices()
pdf(pdfFile, height = height, width = width)
grDevices::pdf(pdfFile, height = height, width = width)
}
index = 0
pb <- progress::progress_bar$new(total = length(plots.l))
for (indexAll in 1:length(plots.l)) {
for (indexAll in seq_len(length(plots.l))) {
pb$tick()
index = index + 1
......@@ -76,7 +76,7 @@
}
if (!checkmate::testNull(pdfFile))
dev.off()
grDevices::dev.off()
futile.logger::flog.info(paste0("Finished writing plots to file ", pdfFile))
}
......@@ -181,8 +181,8 @@
.clearOpenDevices <- function() {
while (length(dev.list()) > 0) {
dev.off()
while (length(grDevices::dev.list()) > 0) {
grDevices::dev.off()
}
}
......@@ -250,13 +250,11 @@
}
}
#' @importFrom BiocParallel multicoreWorkers MulticoreParam
#' @import checkmate
.initBiocParallel <- function(nWorkers, verbose = FALSE) {
#.checkAndLoadPackages(c("BiocParallel"), verbose = verbose)
checkmate::assertInt(nWorkers, lower = 1)
checkmate::assertInt(multicoreWorkers())
checkmate::assertInt(BiocParallel::multicoreWorkers())
if (nWorkers > BiocParallel::multicoreWorkers()) {
warning("Requested ", nWorkers, " CPUs, but only ", BiocParallel::multicoreWorkers(), " are available and can be used.")
......@@ -267,8 +265,7 @@
}
#' @importFrom BiocParallel bplapply bpok multicoreWorkers
#' @import checkmate
.execInParallelGen <- function(nCores, returnAsList = TRUE, listNames = NULL, iteration, abortIfErrorParallel = TRUE, verbose = TRUE, functionName, ...) {
#.checkAndLoadPackages(c("BiocParallel"), verbose = verbose)
......@@ -278,13 +275,13 @@
checkmate::assertFlag(returnAsList)
checkmate::assertFunction(functionName)
checkmate::assertVector(iteration, any.missing = FALSE, min.len = 1)
checkmate::assert(checkmate::checkNull(listNames), checkCharacter(listNames, len = length(iteration)))
checkmate::assert(checkmate::checkNull(listNames), checkmate::checkCharacter(listNames, len = length(iteration)))
res.l = list()
maxCores = tryCatch(
{
out = multicoreWorkers() / 2
out = BiocParallel::multicoreWorkers() / 2
},
error=function() {
......@@ -304,7 +301,7 @@
if (nCores > 1) {
res.l = tryCatch( {
bplapply(iteration, functionName, ..., BPPARAM = .initBiocParallel(nCores))
BiocParallel::bplapply(iteration, functionName, ..., BPPARAM = .initBiocParallel(nCores))
}#, error = function(e) {
# warning("An error occured while executing the function with multiple CPUs. Trying again using only only one CPU...")
......@@ -312,13 +309,13 @@
# }
)
failedTasks = which(!bpok(res.l))
failedTasks = which(!BiocParallel::bpok(res.l))
if (length(failedTasks) > 0) {
warning("At least one task failed while executing in parallel, attempting to rerun those that failed: ",res.l[[failedTasks[1]]])
if (abortIfErrorParallel) stop()
res.l = tryCatch( {
bplapply(iteration, functionName, ..., BPPARAM = .initBiocParallel(nCores), BPREDO = res.l)
BiocParallel::bplapply(iteration, functionName, ..., BPPARAM = .initBiocParallel(nCores), BPREDO = res.l)
}, error = function(e) {
warning("Once again, an error occured while executing the function with multiple CPUs. Trying again using only only one CPU...")
......@@ -372,9 +369,9 @@
hashesStrError = paste0(paste0(rep("#", nchar(lastPartError) - 1), collapse = ""), "\n")
messageError = paste0(objectPart, checkResult, "\n\n", hashesStrError, lastPartError, hashesStrError)
lastPartWarning = "# This warning may or may not be ignored. Carefully check the diagnostic files and downstream results. #\n"
hashesStrWarning = paste0(paste0(rep("#", nchar(lastPartWarning) - 1), collapse = ""), "\n")
messageWarning = paste0(objectPart, checkResult, "\n\n", hashesStrWarning, lastPartWarning, hashesStrWarning)
lastPartWarning = ". \nThis warning may or may not be ignored. Carefully check its significance and whether it may affect the results."
#hashesStrWarning = paste0(paste0(rep("#", nchar(lastPartWarning) - 1), collapse = ""), "\n")
messageWarning = paste0(objectPart, checkResult, lastPartWarning) # , hashesStrWarning)
......@@ -394,7 +391,7 @@
invalidEnd = which(annotation$end > seqlengths[as.character(annotation$chr)] | annotation$end