Pipe SamToFastq, BWA mem not working
1
0
Entering edit mode
4.6 years ago
abbysue ▴ 10

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 :(

Gatk SAM BWA Pipe • 3.4k views
ADD COMMENT
3
Entering edit mode

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

give samtools collate | samtools fastq a try.

ADD REPLY
0
Entering edit mode

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)

ADD REPLY
0
Entering edit mode

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 :)

ADD REPLY
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode

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.

ADD REPLY
1
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

I had the same question!!

ADD REPLY
1
Entering edit mode

Are you on Mac?

ADD REPLY
0
Entering edit mode

Nope (this 20 character min... damn)

ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode
4.3 years ago

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.

ADD COMMENT

Login before adding your answer.

Traffic: 2340 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6