gsnap not producing sam files with header

Question

[Ada]

How can I get gsnap to produce sam files with header? I am having difficulty with current dataset.

[GenoMax]

To get help you need to provide command line you are using. current dataset conveys no information either.

[Ada]

ref=(/dir_001/*.fna.gz)

for f in "${ref[@]}";

do

gmap_build -d genome -D /dir_001 $f

gsnap -D /dir_001 -d genome /R1_001.fastq  /R2_001.fastq -A sam >> "$f.sam"

samtools view -S -b "$f.sam" > "$f.bam";

samtools sort "$f.bam" -o "$f.sorted.bam";

samtools flagstat "$f.sorted.bam" >>  "$f.txt"

done

[Ada]

The code breaks at samtools flagstat, stating header cannot be found.

[GenoMax]

Looks like you forgot to include header when converting SAM to BAM using samtools view.

-h    Include the header in the output.

[Ada]

I made the correction and received this:

[main_samview] fail to read the header from test1.sam file

So it looks like the error is coming from the mapping to output .SAM file command.

The correction of -h to the SAM to BAM file conversion states that their is also failure for reading header in .BAM file.

What can I do to rectify this?

[GenoMax]

What do you get when you try to look at header of one of the files?

samtools -H onefile.sam

Or is there no header at all?

[Ada]

There is no header at all

[GenoMax]

Are you sure the logic of your loop is sound? It looks like you are building a new reference (looping through files) and then trying to align the same fastq files to that index with iteration of the loop.

If you simply try

gmap_build -d genome -D /dir_001 one_fna.fz
gsnap -D /dir_001 -d genome /R1_001.fastq  /R2_001.fastq -A sam > one_fna.sam

does it work properly? If it does then you need to debug your loop. Not sure why you are using >> if you are creating a new SAM files each time.

[Ada]

I am going to debug and follow up

[Ada]

My files were in .fna.gz and that seemed to be the problem in producing SAM file with header.

Converted to .fna and works fine now.