Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • cfrias
    replied
    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

    Leave a comment:


  • swbarnes2
    replied
    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.

    Leave a comment:


  • vaibhavvsk
    replied
    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.

    Leave a comment:


  • Bgansw
    replied
    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

    Leave a comment:


  • Alex Renwick
    replied
    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.

    Leave a comment:


  • lukas1848
    replied
    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.

    Leave a comment:


  • Alex Renwick
    replied
    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

    Leave a comment:


  • maubp
    replied
    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.

    Leave a comment:


  • vaibhavvsk
    replied
    Thanks marcowanger & ulz_peter for your guidelines .It can be easily done by shell,Perl or Java.

    Leave a comment:


  • ulz_peter
    replied
    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!!

    Leave a comment:


  • marcowanger
    replied
    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

    Leave a comment:


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

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-10-2024, 07:30 AM
0 responses
23 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-03-2024, 09:45 AM
0 responses
198 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-03-2024, 08:54 AM
0 responses
209 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-02-2024, 03:00 PM
0 responses
191 views
0 likes
Last Post seqadmin  
Working...
X