Entering edit mode
I am using Diffbind for an ATAC-Seq analysis. My peak caller is MACS2, and here is my sample sheet:
I run Diffbind with the following codes, but it crashed every time on dba.count . it can finished Computing summits... Recentering peaks... Reads will be counted as Paired-end. But have this screen in the end. Does anyone have any idea about this? Thank you so much!
library(BiocParallel)
register(SerialParam())
library(DiffBind)
data <- dba(sampleSheet="sampleSheet.csv")
data <- dba.blacklist(data, blacklist=DBA_BLACKLIST_GRCH38, greylist=FALSE)
data <- dba.count(data)
Error screen:
NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, NA_character_)), elementMetadata = new("DFrame", rownames = NULL, nrows = 117227L, elementType = "ANY", elementMetadata = NULL, metadata = list(), listData = list()), elementType = "ANY", metadata = list()), mode = function (features, reads, ignore.strand = FALSE, inter.feature = TRUE) { features <- .removeSharedRegions(features, ignore.strand = ignore.strand) Union(features, reads, ignore.strand = ignore.strand, inter.feature = inter.feature) }, ignore.strand = TRUE, inter.feature = TRUE, param = new("ScanBamParam", flag = c(keep0 = 4095L, keep1 = 4095L), simpleCigar = FALSE, reverseComplement = FALSE, tag = character(0), tagFilter = list(), what = character(0), which = new("CompressedIRangesList", unlistData = new("IRanges", start = integer(0), width = integer(0), NAMES = NULL, elementType = "ANY", elementMetadata = NULL, metadata = list()), elementType = "IRanges", elementMetadata = NULL, metadata = list(), partitioning = new("PartitioningByEnd", end = integer(0), NAMES = character(0), elementType = "ANY", elementMetadata = NULL, metadata = list())), mapqFilter = 15L), preprocess.reads = NULL), OPTIONS = list( log = FALSE, stop.on.error = TRUE, as.error = TRUE, timeout = NA_integer_, force.GC = FALSE, globalOptions = NULL), BPRNGSEED = c(10407L, 897505460L, 209072368L, -936942767L, -474632773L, 911815027L, 1093500380L), GLOBALS = list(), PACKAGES = character(0))
29: do.call(msg$data$fun, msg$data$args)
30: doTryCatch(return(expr), name, parentenv, handler)
31: tryCatchOne(expr, names, parentenv, handlers[[1L]])
32: tryCatchList(expr, classes, parentenv, handlers)
33: tryCatch({ do.call(msg$data$fun, msg$data$args)}, error = function(e) { list(.error_worker_comm(e, "worker evaluation failed"))})
34: .bpworker_EXEC(msg, bplog(backend$BPPARAM))
35: .recv_any(manager$backend)
36: .recv_any(manager$backend)
37: .manager_recv(manager)
38: .manager_recv(manager)
39: .collect_result(manager, reducer, progress, BPPARAM)
40: .bploop_impl(ITER = ITER, FUN = FUN, ARGS = ARGS, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS, BPREDO = BPREDO, reducer = reducer, progress.length = length(redo_index))
41: bploop.lapply(manager, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS, ...)
42: bploop(manager, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS, ...)
43: .bpinit(manager = manager, X = X, FUN = FUN, ARGS = ARGS, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS, BPREDO = BPREDO)
44: bplapply(X, FUN, ..., BPREDO = BPREDO, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS)
45: bplapply(X, FUN, ..., BPREDO = BPREDO, BPPARAM = BPPARAM, BPOPTIONS = BPOPTIONS)
46: bplapply(setNames(seq_along(reads), names(reads)), function(i, FUN, reads, features, mode, ignore.strand, inter.feature, param, preprocess.reads, ...) { bf <- reads[[i]] .countWithYieldSize(FUN, features, bf, mode, ignore.strand, inter.feature, param, preprocess.reads, ...)}, FUN, reads, features, mode = match.fun(mode), ignore.strand = ignore.strand, inter.feature = inter.feature, param = param, preprocess.reads = preprocess.reads, ...)
47: bplapply(setNames(seq_along(reads), names(reads)), function(i, FUN, reads, features, mode, ignore.strand, inter.feature, param, preprocess.reads, ...) { bf <- reads[[i]] .countWithYieldSize(FUN, features, bf, mode, ignore.strand, inter.feature, param, preprocess.reads, ...)}, FUN, reads, features, mode = match.fun(mode), ignore.strand = ignore.strand, inter.feature = inter.feature, param = param, preprocess.reads = preprocess.reads, ...)
48: .dispatchBamFiles(features, reads, mode, ignore.strand, inter.feature = inter.feature, singleEnd = singleEnd, fragments = fragments, param = param, preprocess.reads = preprocess.reads, ...)
49: .local(features, reads, mode, ignore.strand, ...)
50: GenomicAlignments::summarizeOverlaps(features = intervals, reads = bfl, ignore.strand = TRUE, singleEnd = singleEnd, mode = mode, inter.feature = interfeature, fragments = fragments, param = params)
51: GenomicAlignments::summarizeOverlaps(features = intervals, reads = bfl, ignore.strand = TRUE, singleEnd = singleEnd, mode = mode, inter.feature = interfeature, fragments = fragments, param = params)
52: assay(GenomicAlignments::summarizeOverlaps(features = intervals, reads = bfl, ignore.strand = TRUE, singleEnd = singleEnd, mode = mode, inter.feature = interfeature, fragments = fragments, param = params))
53: pv.getCountsLowMem(bamfile, intervals, bWithoutDupes, mode, yieldSize, singleEnd, fragments, scanbamparam, minCount = minCount, minQC = minMappingQuality, interfeature = interfeature)
54: pv.getCounts(bamfile = countrec$bamfile, intervals = intervals, insertLength = countrec$insert, bWithoutDupes = bWithoutDupes, bLowMem = bLowMem, yieldSize = yieldSize, mode = mode, singleEnd = singleEnd, scanbamparam = scanbamparam, fileType = fileType, summits = summits, fragments = fragments, minMappingQuality = minMappingQuality, minCount = minCount, interfeature = interfeature)
55: FUN(X[[i]], ...)
56: lapply(X = S, FUN = FUN, ...)
57: doTryCatch(return(expr), name, parentenv, handler)
58: tryCatchOne(expr, names, parentenv, handlers[[1L]])
59: tryCatchList(expr, classes, parentenv, handlers)
60: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcall <- deparse(call, nlines = 1L) prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefix, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))})
61: try(lapply(X = S, FUN = FUN, ...), silent = TRUE)
62: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE))
63: FUN(X[[i]], ...)
64: lapply(seq_len(cores), inner.do)
65: mclapply(arglist, fn, ..., mc.preschedule = TRUE, mc.allow.recursive = TRUE)
66: config$lapplyFun(config, params, arglist, fn, ...)
67: dba.parallel.lapply(pv$config, params, todorecs, pv.do_getCounts, bed, bWithoutDupes = bWithoutDupes, bLowMem, yieldSize, mode, singleEnd, scanbamparam, readFormat, summits, fragments, minMappingQuality, minCount, interfeature)
68: pv.counts(res, peaks = res$merged, defaultScore = defaultScore, bLog = bLog, insertLength = insertLength, bOnlyCounts = TRUE, bCalledMasks = TRUE, minMaxval = minMaxval, filterFun = filterFun, bParallel = bParallel, bWithoutDupes = bWithoutDupes, bScaleControl = bScaleControl, bSignal2Noise = bSignal2Noise, bLowMem = saveLowMem, readFormat = readFormat, summits = 0, bRecentered = TRUE, minMappingQuality = minMappingQuality, maxGap = maxGap)
69: pv.counts(DBA, peaks = peaks, minOverlap = minOverlap, defaultScore = score, bLog = bLog, insertLength = fragmentSize, bOnlyCounts = TRUE, bCalledMasks = TRUE, minMaxval = filter, bParallel = bParallel, bUseLast = bUseLast, bWithoutDupes = bRemoveDuplicates, bScaleControl = bScaleControl, filterFun = filterFun, bLowMem = bUseSummarizeOverlaps, readFormat = readFormat, summits = summits, minMappingQuality = mapQCth, minCount = minCount, bSubControl = bSubControl, maxGap = maxGap)
70: dba.count(data)
Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:
My system:
sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.3.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] DiffBind_3.6.1 SummarizedExperiment_1.26.1
[3] Biobase_2.56.0 MatrixGenerics_1.8.0
[5] matrixStats_0.62.0 GenomicRanges_1.48.0
[7] GenomeInfoDb_1.32.2 IRanges_2.30.0
[9] S4Vectors_0.34.0 BiocGenerics_0.42.0
[11] BiocParallel_1.30.3
loaded via a namespace (and not attached):
[1] bitops_1.0-7 RColorBrewer_1.1-3 numDeriv_2016.8-1.1
[4] tools_4.2.0 utf8_1.2.2 R6_2.5.1
[7] irlba_2.3.5 KernSmooth_2.23-20 DBI_1.1.3
[10] colorspace_2.0-3 apeglm_1.18.0 tidyselect_1.1.2
[13] compiler_4.2.0 cli_3.3.0 DelayedArray_0.22.0
[16] rtracklayer_1.56.0 caTools_1.18.2 scales_1.2.0
[19] SQUAREM_2021.1 mvtnorm_1.1-3 mixsqp_0.3-43
[22] stringr_1.4.0 digest_0.6.29 Rsamtools_2.12.0
[25] XVector_0.36.0 jpeg_0.1-9 pkgconfig_2.0.3
[28] htmltools_0.5.2 fastmap_1.1.0 invgamma_1.1
[31] bbmle_1.0.25 limma_3.52.2 BSgenome_1.64.0
[34] htmlwidgets_1.5.4 rlang_1.0.2 BiocIO_1.6.0
[37] generics_0.1.2 hwriter_1.3.2.1 gtools_3.9.2.2
[40] dplyr_1.0.9 RCurl_1.98-1.7 magrittr_2.0.3
[43] GenomeInfoDbData_1.2.8 Matrix_1.4-1 Rcpp_1.0.8.3
[46] munsell_0.5.0 fansi_1.0.3 lifecycle_1.0.1
[49] stringi_1.7.6 yaml_2.3.5 MASS_7.3-57
[52] zlibbioc_1.42.0 gplots_3.1.3 plyr_1.8.7
[55] grid_4.2.0 parallel_4.2.0 ggrepel_0.9.1
[58] bdsmatrix_1.3-6 crayon_1.5.1 lattice_0.20-45
[61] Biostrings_2.64.0 locfit_1.5-9.5 pillar_1.7.0
[64] rjson_0.2.21 systemPipeR_2.2.2 codetools_0.2-18
[67] XML_3.99-0.10 glue_1.6.2 ShortRead_1.54.0
[70] GreyListChIP_1.28.1 latticeExtra_0.6-29 png_0.1-7
[73] vctrs_0.4.1 gtable_0.3.0 purrr_0.3.4
[76] amap_0.8-18 assertthat_0.2.1 ashr_2.2-54
[79] ggplot2_3.3.6 emdbook_1.3.12 restfulr_0.0.15
[82] coda_0.19-4 truncnorm_1.0-8 tibble_3.1.7
[85] GenomicAlignments_1.32.0 ellipsis_0.3.2