Pipe SamToFastq, BWA mem not working

Question

[abbysue]

I'm going through the GATK pipeline to realign bam files. I'm using this command line:

gatk SamToFastq -I markedadapters.bam -F /dev/stdout --CLIPPING_ATTRIBUTE XT --CLIPPING_ACTION 2 --INTERLEAVE true -NON_PF true -TMP_DIR= path/hello | bwa mem -M -t 12 -p /mydata/Homo_sapiens_assembly38.fasta /dev/stdin | gatk MergeBamAlignment --ALIGNED_BAM /dev/stdin --UNMAPPED_BAM unmap.bam -O piped.bam -R /mydata/Homo_sapiens_assembly38.fasta --CREATE_INDEX true --ADD_MATE_CIGAR true --CLIP_ADAPTERS false --CLIP_OVERLAPPING_READS true --INCLUDE_SECONDARY_ALIGNMENTS true --MAX_INSERTIONS_OR_DELETIONS -1 --PRIMARY_ALIGNMENT_STRATEGY MostDistant --ATTRIBUTES_TO_RETAIN XS -TMP_DIR= path/hello

I keep getting this error:

htsjdk.samtools.SAMException: Error in writing fastq file /dev/stdout

I'm about to give up and just run each command separately, but these files are ~150G :(

[Pierre Lindenbaum]

I'm about to give up and just run each command separately,

give samtools collate | samtools fastq a try.

[aalith]

I've read that the behavior of Picard in regressing to a fastq file is preferable (forum post on Biostars, have to find it real quick)

[ATpoint]

I do not see why this would be the case. I am using the above command (and never used Picard for it) all the time, never fails :)

[aalith]

Could abbysue's problem be due to the fact that it's running inside docker and not on the local machine? the piped command listed should work

[WouterDeCoster]

Could you try running the first part before the pipe and see if that works?

gatk SamToFastq -I markedadapters.bam -F /dev/stdout --CLIPPING_ATTRIBUTE XT --CLIPPING_ACTION 2 --INTERLEAVE true -NON_PF true -TMP_DIR= path/hello > mytest.fastq

If one of the steps later fails then the pipe will break. Make sure all components work before making a large piped command.

[abbysue]

I ended up running

gatk SamToFastq -I markedadapters.bam -F f797.fastq --CLIPPING_ATTRIBUTE XT --CLIPPING_ACTION 2 --INTERLEAVE true -NON_PF true -TMP_DIR= path/hello

which of course worked. I tried your command and it worked!! So should I include the " > file.fastq" blurb to get my piped commands to work?

[aalith]

I had the same question!!

[ATpoint]

Are you on Mac?

[abbysue]

Nope (this 20 character min... damn)

[ATpoint]

Ok, on Mac I sometimes had issues with non-writable stdout. As for the character limit, just add whitespaces ;-)

Answer 1

[Marco Pessoa]

What fixed this error for me was to, in addition to indexing the reference with bwa, also indexing the reference using samtools faidx, and creating a sequence dictionary using Picard's CreateSequenceDictionary as a final step.

This (old) tutorial at gatk's forum shows details.