Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • How to extract mapped and unmapped raw reads from bwa's sam file ?

    I have aligned paired end with reference sequence using BWA and got sam output file.I tried picard tools samToFastq tool to convert sam into paired end fastq files but it's for just mapped reads.
    I want to know how to extract unmapped raw reads in separate raw fastq files from sam file.
    Vaibhav Kulkarni

  • #2
    first, identify the unmapped reads from the BAM file using the FLAG field (0x4). Record the name of the read Format file see: http://samtools.sourceforge.net/SAM1.pdf

    Then, read your fastq file, extract the reads matched the names previously recorded
    Marco

    Comment


    • #3
      If you're familiar with scripting, that would be a relatively easy task:
      Unmapped reads have an asterisk (*) in the third column, so you filter out the other ones and then extract the read name (first column) and write it into the output file after a @ sign. Then take the sequence (column 10) write it into the next line; then write the name of the read again after a + sign (it might be enough to just write a + sign). And write the qualities (column 11)in the last row.

      Have fun programming!!

      Comment


      • #4
        Thanks marcowanger & ulz_peter for your guidelines .It can be easily done by shell,Perl or Java.
        Vaibhav Kulkarni

        Comment


        • #5
          Originally posted by ulz_peter View Post
          If you're familiar with scripting, that would be a relatively easy task:
          Unmapped reads have an asterisk (*) in the third column, so you filter out the other ones and then extract the read name (first column) and write it into the output file after a @ sign. Then take the sequence (column 10) write it into the next line; then write the name of the read again after a + sign (it might be enough to just write a + sign). And write the qualities (column 11)in the last row.

          Have fun programming!!
          That may not do exactly what is required. Some unmapped reads don't have an asterisk in the third column (reference name), specifically for paired ends where only one read maps, the other read it given dummy positions to match (for sorting purposes), but the FLAG says it is unmapped.

          i.e. You can only trust the FLAG. One easy way to do this is with the -f and -F arguments to samtools view.

          Comment


          • #6
            Originally posted by maubp View Post
            ...

            i.e. You can only trust the FLAG. One easy way to do this is with the -f and -F arguments to samtools view.
            Yep, that's the way to do it.

            Position 4 indicates an unmapped read. Hence,
            Code:
            samtools view -F4 sample.bam > sample.mapped.sam
            samtools view -f4 sample.bam > sample.unmapped.sam

            Comment


            • #7
              you could also use a simple bash one liner to do that:

              grep ^...............4 aln.sam -v > mapped.aln.sam

              its not pretty, but it works.

              Comment


              • #8
                Originally posted by lukas1848 View Post
                you could also use a simple bash one liner to do that:

                grep ^...............4 aln.sam -v > mapped.aln.sam

                its not pretty, but it works.
                Actually, that won't really work. The flag is a bit field, so the values are additive. If the flag value is 4, then the "unmapped" flag is set and the other ten flags are unset. If a read is unmapped and paired, for example, the flag value would be 4 + 1 = 5.

                Also, there's the less subtle point that your grep expression requires the flag to be exactly 16 characters from the beginning of the line. That won't always be the case.

                Comment


                • #9
                  extract unmapped reads

                  Originally posted by Alex Renwick View Post
                  Yep, that's the way to do it.

                  Position 4 indicates an unmapped read. Hence,
                  Code:
                  samtools view -F4 sample.bam > sample.mapped.sam
                  samtools view -f4 sample.bam > sample.unmapped.sam
                  Hello Alex,

                  I wanted to confirm if I'm interpreting my results correctly.I'm trying to align my sample illumina bacterial genome to an already published reference bacterial genome. The aim is to find out what genes are present in my sample genome, but missing in the reference. The plan of action was to run bwasw, and then pull out the sequence in the sample genome that do not align to the reference, then look at what these sequences are.

                  I aligned my query genome to the ref using bwasw.

                  The output was in sam format. I converted that to bam format using

                  samtools view -bS aln.sam > aln.bam

                  Then I ran the command lines quoted above.

                  When I headed my mapped.sam file, it appears to have worked, however, my unmapped.sam file is empty,the size is zero.
                  This is consistent for two of my sample genomes. It is highly likely that the two sample genomes have no extra genes, when compared to the reference, thus there are no unmapped reads, compared to the reference.

                  However, is there any other way to check this?

                  I'd be very grateful for any help/advice on this.

                  Thanx

                  bgansw

                  Comment


                  • #10
                    Hello Alex,

                    I have also same kind of experience as with Bgansw

                    samtools view -F4 sample.bam > sample.mapped.sam

                    samtools view -f4 sample.bam > sample.unmapped.sam

                    unmapped file is empty.
                    Vaibhav Kulkarni

                    Comment


                    • #11
                      Poking around, I found this thread:

                      Discussion of next-gen sequencing related bioinformatics: resources, algorithms, open source efforts, etc


                      Where someone claims that bwasw doesn't return unmapped reads. You can probably confirm this, by counting reads in your fastq, and reads in the .bam

                      So, you can process the .sam file and the fastq file with a script, to get the read that are in the fastq, but not in the .sam. Or use an aligner that will output the unaligned reads. Bowtie and bwa aln of course for Illumina reads, but you must have longer reads since you aren't using bwa aln, so I'm not sure what's appropriate for those.

                      Comment


                      • #12
                        How to extract multi-mapped reads from bwa's sam file?

                        Hi, I want to extract multi-mapped reads from sam file generated from bwa (bwsw).
                        I am working with 454 data and I have used bwa 0.5.9.

                        The problem is that the flag to indicate that the read is a multiple hit not appear
                        (256 or NH:i: ),
                        So, I can not extract the reads that has multiple hits

                        so I have two question?
                        ## I'ts fine is all flags I can see in the sam file are 0, 4, 16?
                        ## How can I extract the reads that have multiple Hits?

                        e.g: I have two reads mapping in the same contig

                        @SQ SN:Contig1276 LN:500
                        HKU61D301B85SS_nem_043 0 Contig1276 172 27 8S90M1I8M2I18M428S * 0 0 ATCTGTGGAAAAATCATCACCACTGGGGAGTGGAGGAAAGGTAAGATTGCTGAATGGAAGACGAGGTCTCTGCCATTGTCCGATGGCCACATCACCTGCGAAGAATATCAGTTTGAATGTAAAATAGTTTATTGTAAACCTCTGGTTCTCGGATTTTCTTGGTTTCCCTACCATATTAAATGATGGGAAACCAAGAAAATACCAAGAAACGAAGGTTTACAACAAACTGCTTGCAAAAGACGATTTCAGGGTTTTCTGTTTTGCCTTGCAATCACCCCTGTAAGTTCTCATTTATATTTGAATGCTTCACATATTCTTTGTGCATTTTTGTTATATTTTGTTAAATGAGCATGTTTTTTTCCCCCATGCAGAATATGACCAGCAATAACGTTATATGGAATGAGAGTAAGGAGAAAGCAGAAGGAAATCAGTTATCGTCATGGAATGTGGACATTTTTGTGCAAGTTGTGCGAGAAAATGGTGAGATTTTTCACATTATTATTCAGCTTATTTATTAGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGGTGTGGT IIIIIIDD62222DFIIEA??E@A10//42511200115;5400005:;;;::78444457=;<@;888<<<=ACCECDBDEFDDFF@>55<<@330001426?@????79799B@77333990/3----3.44335577:?;;??@<<7999@;;;?@@@>>>:11114/77@@@@A@666??AAB>A3333389>8886/------,1181359?BDDDFFEFIIFFEDDDDFB;;;;;EEEEFB@@@@@?<33325535///58;<;<551225558;;:9;788/- AS:i:76 XS:i:0 XF:i:3 XE:i:1 XN:i:0
                        HKU61D301B85SS_nem_043 0 Contig1276 400 19 204S30M1D46M275S * 0 0 ATCTGTGGAAAAATCATCACCACTGGGGAGTGGAGGAAAGGTAAGATTGCTGAATGGAAGACGAGGTCTCTGCCATTGTCCGATGGCCACATCACCTGCGAAGAATATCAGTTTGAATGTAAAATAGTTTATTGTAAACCTCTGGTTCTCGGATTTTCTTGGTTTCCCTACCATATTAAATGATGGGAAACCAAGAAAATACCAAGAAACGAAGGTTTACAACAAACTGCTTGCAAAAGACGATTTCAGGGTTTTCTGTTTTGCCTTGCAATCACCCCTGTAAGTTCTCATTTATATTTGAATGCTTCACATATTCTTTGTGCATTTTTGTTATATTTTGTTAAATGAGCATGTTTTTTTCCCCCATGCAGAATATGACCAGCAATAACGTTATATGGAATGAGAGTAAGGAGAAAGCAGAAGGAAATCAGTTATCGTCATGGAATGTGGACATTTTTGTGCAAGTTGTGCGAGAAAATGGTGAGATTTTTCACATTATTATTCAGCTTATTTATTAGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGGTGTGGT IIIIIIDD62222DFIIEA??E@A10//42511200115;5400005:;;;::78444457=;<@;888<<<=ACCECDBDEFDDFF@>55<<@330001426?@????79799B@77333990/3----3.44335577:?;;??@<<7999@;;;?@@@>>>:11114/77@@@@A@666??AAB>A3333389>8886/------,1181359?BDDDFFEFIIFFEDDDDFB;;;;;EEEEFB@@@@@?<33325535///58;<;<551225558;;:9;788/- AS:i:57 XS:i:0 XF:i:3 XE:i:1 XN:i:0

                        I think that I can filt it using a Perl script that checking the MAPQ(5 column)?
                        (In that case the higher is 27)
                        It is correct?

                        Thanks in advance,

                        Cris

                        Comment

                        Latest Articles

                        Collapse

                        • seqadmin
                          Exploring the Dynamics of the Tumor Microenvironment
                          by seqadmin




                          The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                          07-08-2024, 03:19 PM
                        • seqadmin
                          Exploring Human Diversity Through Large-Scale Omics
                          by seqadmin


                          In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                          06-25-2024, 06:43 AM

                        ad_right_rmr

                        Collapse

                        News

                        Collapse

                        Topics Statistics Last Post
                        Started by seqadmin, 07-16-2024, 05:49 AM
                        0 responses
                        28 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-15-2024, 06:53 AM
                        0 responses
                        34 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-10-2024, 07:30 AM
                        0 responses
                        41 views
                        0 likes
                        Last Post seqadmin  
                        Started by seqadmin, 07-03-2024, 09:45 AM
                        0 responses
                        205 views
                        0 likes
                        Last Post seqadmin  
                        Working...
                        X