Rnaseq: DupRadar broken

[apeltzer]

Tested v1.4 and it seems that DupRadar is broken:

~/D/t/80791ee1db8a958074d47827688b86> ls -l
total 1508544
-rw-r--r-- 1 alex users 788007786 Apr  9 16:00 AS-225023-LR-34148Aligned.sortedByCoord.out.markDups.bam
-rwxr-xr-x 1 alex users      5285 Apr  9 16:07 dupRadar.r*
-rw-r----- 1 alex users 756721066 Apr  9 16:00 Mus_musculus.GRCm38.90.gtf
alex@aragorn ~/D/t/80791ee1db8a958074d47827688b86> 
singularity shell ../singularity/scilifelab-ngi-rnaseq-1.4.img 
Singularity: Invoking an interactive shell within container...

Singularity scilifelab-ngi-rnaseq-1.4.img:~/Downloads/testme_locally/80791ee1db8a958074d47827688b86> sh .command.sh
Input bam      (Arg 1): AS-225023-LR-34148Aligned.sortedByCoord.out.markDups.bam
Input gtf      (Arg 2): Mus_musculus.GRCm38.90.gtf
Strandness     (Arg 3): unstranded
paired/single  (Arg 4): single
Nb threads     (Arg 5): 2
R package loc. (Arg 6): 2
Output basename       : AS-225023-LR-34148Aligned.sortedByCoord.out.markDups
Loading required package: dupRadar
Loading required package: parallel

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: .C("R_readSummary_wrapper", as.integer(n), as.character(cmd),     PACKAGE = "Rsubread")
 2: Rsubread::featureCounts(files = bam, annot.ext = gtf, isGTFAnnotationFile = TRUE,     GTF.featureType = "exon", GTF.attrType = "gene_id", nthreads = threads,     isPairedEnd = paired, strandSpecific = stranded, ignoreDup = dup,     countMultiMappingReads = mh, ...)
 3: count(mh = TRUE, dup = FALSE)
 4: eval(expr, pf)
 5: eval(expr, pf)
 6: withVisible(eval(expr, pf))
 7: evalVis(expr)
 8: capture.output(counts <- list(mhdup = count(mh = TRUE, dup = FALSE),     mhnodup = count(mh = TRUE, dup = TRUE), nomhdup = count(mh = FALSE,         dup = FALSE), nomhnodup = count(mh = FALSE, dup = TRUE)))
 9: analyzeDuprates(input_bam, annotation_gtf, stranded, paired_end,     threads)
An irrecoverable exception occurred. R is aborting now ...
Segmentation fault (core dumped)
Singularity scilifelab-ngi-rnaseq-1.4.img:~/Downloads/testme_locally/80791ee1db8a958074d47827688b86>

If I change things in the .command.sh to use only 1 CPU core, things work and finish properly.

[apeltzer]

Temporary workaround would be to restrict calling with "2" cores and instead using just one. Runtime for a regular sample with ~2GB BAM file was just a couple of minutes, so I don't really see the need to have everything paralellized. I'm looking into the bioconda recipe, maybe this is just a packaging issue.

[ewels]

Tested on our cluster and everything works fine...

nextflow run SciLifeLab/NGI-RNAseq -r 1.4 -resume -profile uppmax_devel --project snic2017-7-196 --gtf ./test_data/ngi-rna_test_set/genes.gtf --bed12 ./test_data/ngi-rna_test_set/genes.bed --hisat2_index ./test_data/ngi-rna_test_set/r64/ --aligner hisat2 --singleEnd --reads "./test_data/ngi-rna_test_set/*.fastq.gz" -with-singularity scilifelab-ngi-rnaseq-1.4.img

.command.log

/pica/h1/phil/glob/testdata/rnatest/work/b7/b0de8553b39236182887629af3d085/.command.stub: line 38: ps: command not found
/pica/h1/phil/glob/testdata/rnatest/work/b7/b0de8553b39236182887629af3d085/.command.stub: line 27: ps: command not found
Input bam      (Arg 1): SRR4238379_subsamp.sorted.markDups.bam
Input gtf      (Arg 2): genes.gtf
Strandness     (Arg 3): unstranded
paired/single  (Arg 4): single
Nb threads     (Arg 5): 2
R package loc. (Arg 6): 2
Output basename       : SRR4238379_subsamp.sorted.markDups
Loading required package: dupRadar
Loading required package: parallel
null device
          1
null device
          1
null device
          1
[1] "SRR4238379_subsamp.sorted.markDups"

  Sergi SayolsEmail author, Denise Scherzinger and Holger Klein (2016):
  dupRadar: a Bioconductor package for the assessment of PCR artifacts
  in RNA-Seq data. BMC Bioinformatics, 17:428,
  doi:10.1186/s12859-016-1276-2

A BibTeX entry for LaTeX users is

  @Article{,
    title = {dupRadar: a Bioconductor package for the assessment of PCR artifacts in RNA-Seq data},
    author = {Sergi Sayols and Denise Scherzinger and Holger Klein},
    year = {2016},
    journal = {BMC Bioinformatics},
    doi = {10.1186/s12859-016-1276-2},
    url = {http://dx.doi.org/10.1186/s12859-016-1276-2},
    volume = {17},
    issue = {1},
    pages = {428},
  }

R version 3.4.1 (2017-06-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /opt/conda/envs/ngi-rnaseq/lib/R/lib/libRblas.so
LAPACK: /opt/conda/envs/ngi-rnaseq/lib/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  base

other attached packages:
[1] dupRadar_1.8.0

loaded via a namespace (and not attached):
[1] compiler_3.4.1     tools_3.4.1        KernSmooth_2.23-15 Rsubread_1.28.0
/pica/h1/phil/glob/testdata/rnatest/work/b7/b0de8553b39236182887629af3d085/.command.stub: line 12: ps: command not found

[apeltzer]

I just deleted the work directory, and started using nf-core/RNAseq and the problem persists on our cluster infrastructure (binac profile). I'll stick on it and figure out whats wrong...

[apeltzer]

Same issue and I found the problem:

https://support.bioconductor.org/p/104340/#104960

In the old image we used rsubread 1.28.1 from source/ bioconductor which has this issue already fixed.

Bioconda ships the version 1.28.0 which has the bug mentioned above.
Solution: I will update the conda package (which is a no brainer) and then will rebuild the docker image locally and see if this gets things to work again on our side.

As far as I see, it’s not happening on all kinds of data but it might happen from time to time so we should update this.

[ewels]

Fantastic 👍

[apeltzer]

https://github.com/bioconda/bioconda-recipes/pull/8597

[apeltzer]

Local tests were successful! We could simply update the Dockerimage once and things should work again! DupRadar automatically pulls the latest version of Rsubread, so this pulls now 1.28.1 including the fix.

Thanks also to @fbartusch for helping with this ;-)