This is a test version of Biostars. For the public version, visit https://www.biostars.org.
gsnap not producing sam files with header

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

assembly
• 2,097 views • link
updated by Biostar Biostar • written by Ada • Ada

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

• 0 views • link
written by GenoMax GenoMax 
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
• 0 views • link
updated by GenoMax GenoMax • written by Ada • Ada

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

• 0 views • link
written by Ada • Ada

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

-h    Include the header in the output.
• 0 views • link
written by GenoMax GenoMax

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?

• 0 views • link
written by Ada • Ada

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?

• 0 views • link
written by GenoMax GenoMax

There is no header at all

• 0 views • link
written by Ada • Ada

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.

• 0 views • link
written by GenoMax GenoMax

I am going to debug and follow up

• 0 views • link
written by Ada • 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.

• 0 views • link
written by Ada • Ada

0 answers

No answers yet.

Log in to answer this question.

More posts like this