gsnap produces more reads than fastq

[atma_weapon]

hello, i am running gsnap with rna seq data against a reference genome:

gsnap -m 10 -B 5 -t 8 -A sam -d GENOME seqs.fastq > res.sam

but strangely, the output sam file contains more reads than the original fastq files:

$ cat seqs.fastq | echo $((`wc -l`/4))
3776979

$ cat res.sam | grep -v '^ *@' | wc -l
6009141



and I don't understand what that means....

thank you

[twaddlac]

I'm not familiar with gsnap but it seems like there could be a few factors, the most likely one being your grep command. When I try to find out how many alignments have been reported I use the command:

cut -f1 out.sam | sort | uniq | wc -l

That way you don't count nonspecific alignments. That being said, your SAM file may contain multiple alignments for many reads. There should be a flag to control how many alignments per read are reported but if not there are ways around it. I hop this helps!

[severin]

I'm not familiar with gsnap but it seems like there could be a few factors, the most likely one being your grep command. When I try to find out how many alignments have been reported I use the command:

cut -f1 out.sam | sort | uniq | wc -l

That way you don't count nonspecific alignments. That being said, your SAM file may contain multiple alignments for many reads. There should be a flag to control how many alignments per read are reported but if not there are ways around it. I hop this helps!

Besides mapping to multiple locations, the sam file also has a header that can contain many lines and throw off your count.

[twaddlac]

Besides mapping to multiple locations, the sam file also has a header that can contain many lines and throw off your count.

I forgot to mention, to omit the header lines you should do

grep -v '^@' out.sam | cut -f1 | sort | uniq | wc -l

If you convert your SAM file to BAM you will:
A) reduce the size of the file (BAM < SAM)
B) view the alignment output without the headers

samtools view out.bam

To convert the sam file to bam, just do

samtools view -bS out.sam > out.bam

[Simon Anders]

If a read is aligned to multiple location, it appears multiple times in the SAm file. Sort the file by read name, cut to the first coulmn (read name), pipe it through 'uniq', then count again.