Find unmapped read from sam/bam file

[genelab]

Hi guys

Could anyone tell me how to extract the unmapped reads from sam/bam file and convert the result to fastq format using samtools?

Thx!

[Adrian_H]

samtools view -f 4 yourbamfile.bam will give you unmapped reads

Then pull out the first column of read names (cut -f1 -d" ") and extract those reads from your original fastq files, or make an awk script to reformat the readid, sequence, and quality scores into fastq.

Note that depending on the alignment program that you are using, unmapped reads may or may not be reported in the results. Also, some programs trim off the /1 or /2 of the read ID, if you are working with paired ends). Finally, keep in mind that if you use this to extract other flags, the sequence in the BAM file is only what aligned, and could be the reverse complement of the input. (Shouldn't be an issue for unmapped reads)

[Adamo]

If you've used bwasw then the command line suggested by Adrian won't work.
You'll have to write your own perl or awk script to extract unmapped reads comparing your bam output with your fastq file and rewrite it omitting aligned reads.

[bhootnaath]

can use bam2fastq

http://www.hudsonalpha.org/gsl/software/bam2fastq.php

[csquared]

+1 on the BAM2FASTQ. Great tool...of course I'm biased as it came from my group but it is well documented and fast. Let us know if you have any questions or problems.

[byb121]

Hi,

I used bam2fastq tool to extract unmapped reads, it's really fast and better documented. but I had difficulties to address the cause of the warning message:

$ ./bam2fastq -o s_%#_extracted_reads.txt -f --no-aligned --unaligned --no-filter alignments.bam 
[bam_header_read] EOF marker is absent.
This looks like paired data from lane 1.
Output will be in s_1_1_extracted_reads.txt and s_1_2_extracted_reads.txt
55130926 sequences in the BAM file
8238703 sequences exported
WARNING: 5947209 reads could not be matched to a mate and were not exported

Fastq files contain 1145747 reads each, which means those 5947209 unmapped reads are discarded. But I really would like to have them included in the result. Could you help me out here?

PS: Reads are pair end, ranging from 25 - 78 after quality trimming.

+1 on the BAM2FASTQ. Great tool...of course I'm biased as it came from my group but it is well documented and fast. Let us know if you have any questions or problems.

[vishal.rossi]

samtools view -bh -f 0*4 -o output.file input

[JonB]

What about reads mapping to the reverse strand? Should they be reverse complemented before converting to fastq?

[Brian Bushnell]

What about reads mapping to the reverse strand? Should they be reverse complemented before converting to fastq?

No, sequences and qualities are always the same as the source fastq, regardless of mapped strand.

[JonB]

Good to know, thanks!