Hi, I have tried to run R package- DiffBind in SSH (Linux). I have R script works will with windows R 3.3.3 but failed at SSH sever linux R. It stops at dba.count step. Could anyone help me? Does it mean the bam file is required instead of bed file? Thank you.
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
library(DiffBind)
setwd("/net/wonderland/home/zhaolinz/practice/batch8_macs2")
#list.files()
#list.files("diff")
diff <- dba(sampleSheet="CLA050417.csv",minOverlap=2,config=data.frame(RunParallel=TRUE, reportInit="diff",AnalysisMethod=DBA_DESEQ2,minQCth=5,fragmentSize=300,bCorPlot=FALSE,th=0.1, bUsePval=FALSE))
diff_count <- dba.count(diff, minOverlap=2,score=DBA_SCORE_TMM_READS_FULL, fragmentSize=diff$config$fragmentSize, filter=0, bRemoveDuplicates=FALSE, filterFun=max,bCorPlot=diff$config$bCorPlot,
bUseSummarizeOverlaps=FALSE, readFormat=DBA_READS_DEFAULT,bParallel=diff$config$RunParallel)
.....
<font face="monospace">The error is here </font>
*** caught segfault ***
address (nil), cause 'unknown'
Traceback:
1: .Call("croi_count_reads", bamfile, as.integer(insertLength), as.integer(fileType), as.integer(bufferSize), as.integer(minMappingQual), as.character(intervals[[1]]), as.integer(intervals[[2]]), as.integer(intervals[[3]]), as.integer(icount), as.logical(bWithoutDupes), as.logical(bSummits), counts, summits.vec, heights.vec)
2: cpp_count_reads(bamfile, insertLength, fileType, bufferSize, intervals, bWithoutDupes, summits, minMappingQuality)
3: 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)
4: FUN(X[[i]], ...)
5: lapply(X = S, FUN = FUN, ...)
6: doTryCatch(return(expr), name, parentenv, handler)
7: tryCatchOne(expr, names, parentenv, handlers[[1L]])
8: tryCatchList(expr, classes, parentenv, handlers)
9: 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)[1L] pr
efix <- paste("Error in", dcall, ": ") LONG <- 75L msg <- conditionMessage(e) sm <- strsplit(msg, "\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, conditio
nMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && identical(getOption("show.error.messages"), TRUE)) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))})
10: try(lapply(X = S, FUN = FUN, ...), silent = TRUE)
11: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE))
12: FUN(X[[i]], ...)
13: lapply(seq_len(cores), inner.do)
14: mclapply(arglist, fn, ..., mc.preschedule = TRUE, mc.allow.recursive = TRUE)
15: config$lapplyFun(config, params, arglist, fn, ...)
16: dba.parallel.lapply(pv$config, params, todorecs, pv.do_getCounts, bed, bWithoutDupes = bWithoutDupes, bLowMem, yieldSize, mode, singleEnd, scanbamparam, readFormat, summits, fragments, minMappingQuality)
17: pv.counts(DBA, peaks = peaks, minOverlap = minOverlap, defaultScore = score, bLog = bLog, insertLength = fragmentSize, bOnlyCounts = T, bCalledMasks = TRUE, minMaxval = filter, bParallel = bParallel, bUseLast = bUseLast, bWithoutDupes = bRemoveDuplicates, bScaleControl = bScaleControl, filterFun = filterFun, bLowMem = bUseSummarizeOverlaps, readFormat = readFormat, summits = summits, minMappingQuality = mapQCth)
18: dba.count(diff, minOverlap = 3)

Hi,
Could you post the output of "
sessionInfo()" so we can see what version you're running? There were bugs in earlier versions of DiffBind when using bed files for counts. Just to clarify your comment, are your reads in.bed(or.bed.gz) format instead of.bam? We used to support.bedfor reads as well as peaks, but we decided to remove it (though I can't off-hand recall whether or not I actually did take the code for it). And if your reads are indeed.bed, post the first few lines of one of them.(As a side note, if you tag your post with 'diffbind' we'll get email when you post, so you might get a faster response.)
Your sample sheet looks reasonable. (Normally we'd store all the bam or bed or peak files in subdirectories called
bamorbedorpeak, so the sample sheet would have to include the full (or relative) paths to the files, but you would get a different error if it couldn't find the files at all.)Please post the output of the
sessionInfo()command, and the first few lines of one of your .bed files, so I can see whether the old bug is the problem.(It turns out that we do still support
.bedas a format for reads.)- Gord
It works when DBA was built but error stopped at dba.count(). I had ran the same script using R program in PC and it worked fine. However, short of memory space led me to follow the SSH linux system.
300073_57756 CD4 CLA Positive 0h 9 narrow
300073_57757 CD4 CLA Negative 0h 9 narrow
300073_57758 CD8 CLA Positive 0h 9 narrow
300073_57759 CD8 CLA Negative 0h 9 narrow
*** caught segfault ***
address (nil), cause 'unknown'
Traceback:
1: .Call("croi_count_reads", bamfile, as.integer(insertLength), as.integer(fileType), as.integer(bufferSize), as.integer(minMappingQual), as.character(intervals[[1]]), as.integer(intervals[[2]]), as.integer(intervals[[3]]), as.integer(icount), as.logical(bWithoutDupes), as.logical(bSummits), counts, summits.vec, heights.vec)
There is a command in R called "
sessionInfo()" that shows what versions of R and the currently-loaded packages you are running. For example, on a slightly out-of-date version R, if I execute the command "sessionInfo()" at the R prompt, I get these results:Could you (after running the command "
library(DiffBind)" ) please type the command "sessionInfo()" at the R command prompt and post the results?After you've done that, could you, at the SSH command prompt, i.e. the shell command line, type
head 64502_macs2_peaks.narrowPeakor within an R session, try
system("head 64502_macs2_peaks.narrowPeak")so that I can see the first few lines of one of your read files?
I've asked you for these things twice before; this is the third time. I can't diagnose the problem without this information.
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.5 LTS
Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0
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 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] DiffBind_1.16.3 RSQLite_1.1-2
[3] locfit_1.5-9.1 GenomicAlignments_1.6.3
[5] Rsamtools_1.22.0 Biostrings_2.38.4
[7] XVector_0.10.0 limma_3.26.9
[9] SummarizedExperiment_1.0.2 Biobase_2.30.0
[11] GenomicRanges_1.22.4 GenomeInfoDb_1.6.3
[13] IRanges_2.4.8 S4Vectors_0.8.11
[15] BiocGenerics_0.16.1
loaded via a namespace (and not attached):
[1] Rcpp_0.12.10 lattice_0.20-35 GO.db_3.2.2
[4] gtools_3.5.0 digest_0.6.12 plyr_1.8.4
[7] backports_1.0.5 futile.options_1.0.0 BatchJobs_1.6
[10] ShortRead_1.28.0 ggplot2_2.2.1 gplots_3.0.1
[13] zlibbioc_1.16.0 GenomicFeatures_1.22.13 lazyeval_0.2.0
[16] annotate_1.48.0 gdata_2.17.0 Matrix_1.2-10
[19] checkmate_1.8.2 systemPipeR_1.4.8 GOstats_2.36.0
[22] splines_3.4.0 BiocParallel_1.4.3 stringr_1.2.0
[25] pheatmap_1.0.8 RCurl_1.95-4.8 biomaRt_2.26.1
[28] munsell_0.4.3 sendmailR_1.2-1 compiler_3.4.0
[31] rtracklayer_1.30.4 base64enc_0.1-3 BBmisc_1.11
[34] fail_1.3 tibble_1.3.0 edgeR_3.12.1
[37] XML_3.98-1.7 AnnotationForge_1.12.2 bitops_1.0-6
[40] grid_3.4.0 RBGL_1.46.0 xtable_1.8-2
[43] GSEABase_1.32.0 gtable_0.2.0 DBI_0.6-1
[46] magrittr_1.5 scales_0.4.1 graph_1.48.0
[49] KernSmooth_2.23-15 amap_0.8-14 stringi_1.1.5
[52] hwriter_1.3.2 genefilter_1.52.1 latticeExtra_0.6-28
[55] futile.logger_1.4.3 brew_1.0-6 rjson_0.2.15
[58] lambda.r_1.1.9 RColorBrewer_1.1-2 tools_3.4.0
[61] Category_2.36.0 survival_2.41-3 AnnotationDbi_1.32.3
[64] colorspace_1.3-2 caTools_1.17.1 memoise_1.1.0
> system("head 64502_macs2_peaks.narrowPeak")
1 10340 10634 64502_macs2_peak_1 59 . 4.99269 8.081125.94287 187
1 565190 565464 64502_macs2_peak_2 759 . 8.94809 79.14370
75.99088 155
1 566790 567199 64502_macs2_peak_3 80 . 3.00136 10.24183
8.04778 272
1 569261 569564 64502_macs2_peak_4 937 . 10.13605 97.11760 93.74061 188
1 713661 714645 64502_macs2_peak_5 1063 . 16.07595 109.91768 106.36751 549
1 752597 752887 64502_macs2_peak_6 34 . 3.92283 5.496043.44459 118
1 762453 763263 64502_macs2_peak_7 578 . 14.14343 60.75053 57.81964 410
1 779676 780263 64502_macs2_peak_8 290 . 11.01695 31.59573 29.03273 391
1 793444 793891 64502_macs2_peak_9 77 . 5.70594 9.955587.76720 111
1 794057 794275 64502_macs2_peak_10 59 . 4.99269 8.081125.94287 133
Sorry for the late response. We are updating all the internet system within the department.