We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

SplitNCigarReads fails on unmapped reads in second pass (Illegal argument exception)

Hi,
I'm trying to get RNA-seq data analysis-ready for variant calling according to Best Practices (RNAseq short variant discovery (SNPs + Indels)). SplitNCigarReads is failing and it seems like it might be a bug - I have pored over the previous steps and ValidateSamFile gives the bam file resulting from MarkDuplicates a clean bill of health. Also there have been a couple of similar questions on forums over the years but they were not resolved or closed as 'stale and probably fixed by now', i.e. https://github.com/bcbio/bcbio-nextgen/issues/2354 and https://gatkforums.broadinstitute.org/gatk/discussion/12547/splitncigarreads-fails-on-illegalargumentexception-contig-must-be-non-null-and-not-equal-to-and

Details:
Data - RNA-seq data downloaded from SRA, (SRR5511204); paired-end Illumina, 101bp reads. Alignment done with 2-pass STAR.
Platform: running locally from downloaded GATK version 4.1.4.0 on Ubuntu 16.04
Command (copied from program output):

```
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/kate/Downloads/GATK/gatk-4.1.4.0/gatk-package-4.1.4.0-local.jar SplitNCigarReads -I marked_duplicates2.bam -O splitNcigarred.bam -R /home/kate/Databases/GRCh38/GRCh38.primary_assembly.genome.fa
```
Error (which seems to happen when program gets to unmapped reads in second pass):

```
java.lang.IllegalArgumentException: contig must be non-null and not equal to *, and start must be >= 1
at org.broadinstitute.hellbender.utils.read.SAMRecordToGATKReadAdapter.setMatePosition(SAMRecordToGATKReadAdapter.java:197)
at org.broadinstitute.hellbender.tools.walkers.rnaseq.OverhangFixingManager.setPredictedMateInformation(OverhangFixingManager.java:445)
at org.broadinstitute.hellbender.tools.walkers.rnaseq.SplitNCigarReads.splitNCigarRead(SplitNCigarReads.java:222)
at org.broadinstitute.hellbender.tools.walkers.rnaseq.SplitNCigarReads.secondPassApply(SplitNCigarReads.java:185)
at org.broadinstitute.hellbender.engine.TwoPassReadWalker.lambda$traverseReads$0(TwoPassReadWalker.java:72)
at java.util.stream.ForEachOps$ForEachOp$OfRef.accept(ForEachOps.java:184)
at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
at java.util.stream.ReferencePipeline$2$1.accept(ReferencePipeline.java:175)
at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
at java.util.Iterator.forEachRemaining(Iterator.java:116)
at java.util.Spliterators$IteratorSpliterator.forEachRemaining(Spliterators.java:1801)
at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:481)
at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:471)
at java.util.stream.ForEachOps$ForEachOp.evaluateSequential(ForEachOps.java:151)
at java.util.stream.ForEachOps$ForEachOp$OfRef.evaluateSequential(ForEachOps.java:174)
at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:234)
at java.util.stream.ReferencePipeline.forEach(ReferencePipeline.java:418)
at org.broadinstitute.hellbender.engine.TwoPassReadWalker.traverseReads(TwoPassReadWalker.java:70)
at org.broadinstitute.hellbender.engine.TwoPassReadWalker.traverse(TwoPassReadWalker.java:59)
at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1048)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
at org.broadinstitute.hellbender.Main.main(Main.java:292)
```

Hope you can help because I'd really like to try variant calling!
Thanks :)

Answers

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @kescullator

    Would you please share your input and reference files so we can test this on our end? Please follow this link for information on sharing your data with us: https://software.broadinstitute.org/gatk/guide/article?id=1894

  • kescullatorkescullator AustraliaMember
    Sorry Bhanu,
    I'm not familiar with ftp. I tried to upload my archive via opening an ftp connection but it said '425 Unable to build data connection: Connection timed out'. I also tried with curl
    (curl -T kescullator_gatk.tar.gz ftp://ftp.broadinstitute.org --user gsapubftp:5WvQWSfi)
    which failed after 3 sec with this error:
    curl: (25) Failed FTP upload: 550
    Can you suggest what I'm doing wrong?
    Thanks :)
  • mshandmshand Member, Broadie, Dev ✭✭✭

    Hi @kescullator,

    I just tried to look at the SRA data you mentioned, but I'm not familiar with SRA and I'm having trouble determining what format those files are. Are they aligned or did you do the alignment with STAR yourself? Did you use MergeBamAlignment after STAR?

    It's very possible this is still a bug in SplitNCigarReads, but it would be helpful to get your files through the FTP so we can test with them directly. @bhanuGandham do you know what the issue with the ftp connection is?

  • kescullatorkescullator AustraliaMember
    Hi @mshand,
    Thanks for your comment. I downloaded the files from SRA as fastq files using the SRA toolkit fasterq-dump tool. Then I used cutadapt to get rid of the adapter sequences and aligned the reads with STAR using GRCh38 reference genome. Next, trying to follow GATK best practices, I sorted the alignment sam file by queryname (SortSam), and got a uBAM by running the post-cutadapt fastq files through FastqToSam (with -RG, -SM and -PL options), so I could use MergeBamAlignment. Then MarkDuplicates, and then it broke when I tried SplitNCigarReads.

    It would be easier to make sure we're all looking at the same thing by getting my files to you. If ftp isn't working, can we do something else? The tar.gz is 3.4GB, can I share a dropbox folder with you or share via cloudstor or something? It's entirely possible I'm just doing something wrong with ftp...
    Thanks for your help :)
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @kescullator

    I just tested it out myself and the ftp server is working fine. Take a look at these docs for more info on how to use transfer files using ftp:https://tecadmin.net/download-upload-files-using-ftp-command-line/
    Another option is to use a ftp client.

    For more questions on the use of ftp, I think biostars or seqanswers might be a better forum. Those forums might already have ftp usage related questions answered.

  • kescullatorkescullator AustraliaMember
    Thanks @bhanuGandham; I'm sure that was what I was doing... I just tried again but this time it said 'overwrite permission denied'! So does that mean it actually already transferred? It's called kescullator_gatk.tar.gz.
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @kescullator

    Can you please share the command you are using?

  • kescullatorkescullator AustraliaMember
    On cli:
    connect with
    ftp ftp.broadinstitute.org
    asks for name - I give gsapubftp
    password - 5WvQWSfi
    as per https://software.broadinstitute.org/gatk/documentation/article.php?id=1894
    It connects.
    Then I do:
    put kescullator_gatk.tar.gz

    It fails. I can't see many places for error there. Is it possible my institution is blocking it because ftp is not secure? They're very tight with security especially since Australian National University was hacked!
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @kescullator before the put command, try
    mkdir kescullator
    put kescullator_gatk.tar.gz kescullator/

  • kescullatorkescullator AustraliaMember
    No worries - This is what happened.

    ftp ftp.broadinstitute.org
    Connected to ftp01.broadinstitute.org.
    220 FTP Server ready.
    Name (ftp.broadinstitute.org:kate): gsapubftp
    331 Password required for gsapubftp
    Password:
    230 Anonymous access granted, restrictions apply
    Remote system type is UNIX.
    Using binary mode to transfer files.
    ftp> mkdir kescullator
    257 "/kescullator" - Directory successfully created
    ftp> put kescullator_gatk.tar.gz kescullator/
    local: kescullator_gatk.tar.gz remote: kescullator/
    200 PORT command successful
    550 kescullator/: Overwrite permission denied
    ftp>
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin
    edited December 2019

    @kescullator

    I am not sure why this is happening. These commands seem to be working fine on my end. This question might be better suited for biostars/stackoverflow/seqanswers.

Sign In or Register to comment.