Why do I lose sequences when converting my bam files into fastq?

Question

[DNAngel]

Hi all,

I used bamToFastq to convert my bamfiles into R1 and R2 fastq files. I tried samtools fastq before but found it was erroneous (so many missing sequences). However, I still have an issue with bamTofastq from bedtools. When I add the read counts obtained from grep -c "@" R1.fastq and grep -c "@" R2.fastq, it is always slightly less than the count from the bamfile samtools view -c in.bamfile. Why might this be the case? I haven't found any documentation to suggest that this is a normal thing. The R1 and R2 fastq counts should equal the counts in the bamfile so what am I doing wrong with the conversion??

Thank you.

[Ram]

@ is a valid quality encoding character in FASTQ QUAL lines, so you may want to use the simpler line_count / 4 formula to count number of reads in the FASTQ.

[GenoMax]

Or simply add whatever is next to @ in your grep. e.g. @A00500 (generally sequencer serial). That will only count line 1 from each fastq record.

[DNAngel]

Yea I tried both ways, same thing. Just less reads in my fastq files than in my bamfiles.

Answer 1

[Pierre Lindenbaum]

it is always slightly less than the count from the bamfile samtools view -c in.bamfile. Why might this be the case?

you bam contains supplementary, secondary alignments.

[ATpoint]

https://www.biostars.org/p/9475071/ silly me :-P

wanted to link:

losing reads bam to fastq

[DNAngel]

Is this something that I can safely ignore then? I used sambamba to get my bamfiles because samtools was just too unbelievably slow.

[DNAngel]

Yes you are right. There are secondary alignments present. How does bamToFastq know then to not convert these reads into fastq reads??

Answer 2

[swbarnes2]

Are you absolutely sure that every single read in the bam exists in one and only one line? That you have zero supplemental or secondary reads?

[dsull]

Agreed with this response. I have encountered this in my work; e.g. after STAR alignment, my bam file contains multimapping (e.g. a read that aligns to two locations could appear twice or more in the bam file).

Why don't you try samtools view in.bamfile | cut -f1 | sort -u | wc -l to get a count of the number of unique read names in your bam file?

[DNAngel]

Tried this but it cuts out like half of the sequences. I used sambamba to obtain my bamfiles since it is faster than samtools. For the mapped reads I used sambamba view -f bam -F 'not (unmapped or mate_is_unmapped)' in.bam and then the opposite for unmapped reads (just took the not out). samtools view -c in.bam for example gives 436840 reads, but your command gives 218264.

[DNAngel]

I checked, I do have secondary alignments. So I am now curious to know how bamToFastq works. When I check for my reads and exclude secondary alignments in my bamfile, the number makes sense. It matches the number of R1 and R2 fastq reads. How does bamToFastq ignore secondary alignments during conversion? I just want to be sure it is okay and I can continue with my fastq files confident that no secondary alignments were included.

[swbarnes2]

It's likely ignoring all reads which are flagged as secondary and supplemental. This wouldn't work right if you messed around with your bam and removed primary alignments. If you really want to be sure, grep out and count the # of unique read names.