Announcement

Collapse
No announcement yet.

Find unmapped read from sam/bam file

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Find unmapped read from sam/bam file

    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!

  • #2
    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)
    Last edited by Adrian_H; 07-01-2010, 07:52 AM. Reason: changed filename

    Comment


    • #3
      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.

      Comment


      • #4
        can use bam2fastq

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

        Comment


        • #5
          +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.
          HudsonAlpha Institute for Biotechnology
          http://www.hudsonalpha.org/gsl

          Comment


          • #6
            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:

            Code:
            $ ./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.


            Originally posted by csquared View Post
            +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.

            Comment


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

              Comment


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

                Comment


                • #9
                  Originally posted by JonB View Post
                  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.

                  Comment


                  • #10
                    Good to know, thanks!

                    Comment

                    Working...
                    X