Qualimap RNAseq failed to create plot after BAM file analysis is over

[Haroon Kalam]

Hi ,
I was trying to run an rnaseq module from Qualimap command line and received an error as follow :


BAM file analysis finished
Creating plots
Failed to run rnaseq

java.lang.ArithmeticException: / by zero
    at org.bioinfo.ngs.qc.qualimap.common.TranscriptDataHandler.createJunctionAnalysisPieChart(TranscriptDataHandler.java:474)
    at org.bioinfo.ngs.qc.qualimap.common.TranscriptDataHandler.createPlots(TranscriptDataHandler.java:549)
    at org.bioinfo.ngs.qc.qualimap.process.RNASeqQCAnalysis.createCharts(RNASeqQCAnalysis.java:105)
    at org.bioinfo.ngs.qc.qualimap.process.RNASeqQCAnalysis.createResultReport(RNASeqQCAnalysis.java:92)
    at org.bioinfo.ngs.qc.qualimap.process.RNASeqQCAnalysis.run(RNASeqQCAnalysis.java:69)
    at org.bioinfo.ngs.qc.qualimap.main.RnaSeqQcTool.execute(RnaSeqQcTool.java:135)
    at org.bioinfo.ngs.qc.qualimap.main.NgsSmartTool.run(NgsSmartTool.java:187)
    at org.bioinfo.ngs.qc.qualimap.main.NgsSmartMain.main(NgsSmartMain.java:103)




I am on a 64bit Linux machine and using Qualimap_v2.0 with R version 3.0.2

Any help would be appreciated

[larry.wi...@gmail.com]

Hi Haroon,

Any luck on this problem from the Qualimap developers? I'm having the exact same issue. I've worked through all of the other issues in getting Qualimap rnaseq to work with a human gtf file. I've tried both development versions supplied graciously by Konstantin Okonechnikov in this thread: https://groups.google.com/forum/#!topic/qualimap/M2U9aLVcGj4, as well as the released version. When I say I got QM to work with the gtf file, I mean that it now apparently processes the gtf but fails after the 'creating plots' step, as you show above. It's disappointing, the reports are beautiful.

I've examined my gtf files for all the likely culprits. It would be super useful if the QM people could post a human gtf file that works with their software, or vice-versa.

Larry Wilhelm
Bioinformatics / OHSU

[Haroon Kalam]

Hello Larry,

I couldn't solve this issue but I did a round about by writing small scripts to get statistics from BAM files and then plotted them using ggplot2 package in R.

I believe the error is due to incompatibility between R version and corresponding java call.

Regards

Haroon Kalam

Doctoral Fellow,

Dr. Dhiraj Kumar Lab,

Immunology Research Group,

International Centre For Genetic Engineering And Biotechnology, New Delhi.

--
You received this message because you are subscribed to a topic in the Google Groups "QualiMap" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/qualimap/edvvSaWp-yk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to qualimap+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Konstantin Okonechnikov]

Dear Haroon and Larry,

thanks a lot for the reports!

Actually RNA-seq QC doesn't depend on R, only Counts QC and Multisample BAM QC depend on it. It is some Java code bug.

I would like to take a look on the data to fix this problem.
Could you please run with the small BAM subfile and check if it also fails?

Subsample could be created using samtools:

samtools view -s 0.01 -b file.bam > subsample.bam

Please let me know if Qualimap crashes on the subsample as well, and if it happens please share the subsample file and the GTF. Additionally, you could share the large BAM file by sending link to my e-mail, I will only use it to fix the issues in Quailmap.

Additionally: Larry, thanks a lot for suggestion about GTF files, I will take this into account.

Best regards,

   Konstantin

Best regards,
   Konstantin

--
You received this message because you are subscribed to the Google Groups "QualiMap" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qualimap+u...@googlegroups.com.

[Haroon Kalam]

Dear Konstantin,

After your suggestion I tried RNAseq QC on two different bam files

Case-1. Bam file generated from Tophat with -G option (i.e. with gtf file supplied during alignment)

Case-2. Bam file generated from Tophat without -G option (i.e. without any gtf file)

I can replicate the error only in the case where no gtf file was used during tophat alignment.

GTF file was downloded from " ftp://ftp.ensembl.org/pub/release-78/gtf/homo_sapiens/ " and paired end reads were aligned to GRCh38 using Tophat v2.0.12
I am on a 64bit Linux machine and using Qualimap_v2.0 with java version "1.7.0_75"

Regards

Haroon Kalam

Doctoral Fellow,

Dr. Dhiraj Kumar Lab,

Immunology Research Group,

International Centre For Genetic Engineering And Biotechnology, New Delhi.

[Konstantin Okonechnikov]

Dear Haroon,

thanks a lot for the info! I have access to some human RNA-seq datasets aligned with Tophat, I will try checking it with Qualimap RNAseq QC to detect the bug.

Anyway, it would be great if you could share a subsample from your data where Qualimap fails.

Thanks,

   Konstantin

[larry.wi...@gmail.com]

Hi Konstantin and Haroon,

Thanks for looking into this. I'm not re-creating the same error after subsampling the bam, but a different one. I'm working on fully characterizing the issue on my end and will send a bam and gtf that display the problem as soon as I can (work has me busy).

Larry

[Konstantin Okonechnikov]

Dear Larry,

thanks a lot for your support! Look forward to check your data and fix the problem :)

--

  Konstantin

[Konstantin Okonechnikov]

Dear Haroon,

I' ve detected and fixed the problem with loading the GTF file from Ensembl release-78. The issue was because of semicolon in some gene names:

For example PRAMEF6:
 1    ensembl    gene    12824942    12828663    .    -    .    gene_id "ENSG00000279195"; gene_version "1"; gene_name "PRAMEF6;"; gene_source "ensembl"; gene_biotype "protein_coding";

The build with this fix can be downloaded from here:
https://bitbucket.org/kokonech/qualimap/downloads/qualimap-build-06-02-15.tar.gz

Additionally it's important to allow more RAM to apply RNA-seq QC with this huge GTF file. Here is an example command:

~/qualimap-qualimap/qualimap-build-06-02-15/qualimap rnaseq -bam VCap500.accepted_hits.bam -gtf /data/annotations/Homo_sapiens.GRCh38.78.gtf --java-mem-size=6G

I will be glad to know if there are any other issues.

--

  Konstantin

On Tue, Feb 3, 2015 at 3:48 PM, Haroon Kalam <haroon....@gmail.com> wrote:

[François Lefebvre]

Hello Constantin, I am using build-06-02-15  but a (similar) error still occurs. Log and data to reproduce is here:

https://www.dropbox.com/s/pnc4zyjo3i9gvhi/qualimap_bug_Report.zip

Cheers

Francois L.

[Konstantin Okonechnikov]

Dear Francois,

thanks a lot for the report! I've fixed the issue. The problem was that the BAM file doesn't include any junctions even though this is RNA-seq data. Now if the junctions are not detected the plot with junctions analysis is skipped.

The fixed version is available from here:
https://bitbucket.org/kokonech/qualimap/downloads/qualimap-build-09-02-15.tar.gz

The junctions are usually reported in BAM format using N CIGAR operation (for example as in Tophat). By some reason they are not reported by the tool you used (Butter), or maybe reported somehow different. Do you know these details? Additionally by some reason 5' bias and 3' bias are zeroes, this is also interesting.

Also, let me know if there are any other problems.

--

 Konstantin

[François Lefebvre]

Hi Konstantin, this data is small RNA-seq aligned with butter (=bowtie1) which explains the absence of junctions and bias. Most of the fragments sequenced are full ~21nt miRNAs.

Thanks for Qualimap!

Francois L.

François Lefebvre

[Larry Wilhelm]

Hi Konstanin,

I'm happy to report that I've got qualimap2.0 working, using a bam file created with GSNAP version 2014-08-04, and this GTF file, right out of the box:

ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

My previous issue was due to how the bam file was constructed. I suspected this, your message to Francois regarding the presences of splice junctions and cigar strings confirmed it, and once I found the time to regenerate them with GSNAP, all is well. Works great. 

Thank you!!

Larry Wilhelm OHSU/Bioinformatics

[Konstantin Okonechnikov]

Dear Francois and Larry,

thanks a lot for your quick replies! It is very useful that you reported about problems!

Let me know if there are any other suggestions :)

--

  Konstantin

[정웅재]

Hi, Francois!

I wonder if you could teach me what does aligned with butter means.

I have used qualimap to QC smRNA-seq data of mine, and some samples don't have the 5' 3' bias values. I hope your comment could clear me out why these thing happen in some samples.

Our data is single-end sequencing if this information is needed.

Thanks!

2015년 2월 10일 화요일 오전 3시 55분 11초 UTC+9, François Lefebvre 님의 말: