SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to get uniquely mapped reads from Tophat subeet Bioinformatics 10 11-28-2012 05:56 AM
not uniquely mapped reads unidodo RNA Sequencing 2 04-22-2011 02:07 PM
Uniquely mapped reads with bowtie mapper Bioinformatics 2 11-22-2010 10:44 PM
BWA Uniquely Mapped Reads NF_seq Bioinformatics 0 09-06-2010 03:32 AM
cufflinks and non-uniquely mapped reads clariet Bioinformatics 1 05-08-2010 11:13 AM

Reply
 
Thread Tools
Old 04-24-2012, 04:50 PM   #1
hanleng
Member
 
Location: Houston, TX

Join Date: Mar 2012
Posts: 15
Default Remove reads which are not uniquely mapped

I have already had the BAM file, so how can I remove those reads which are not uniquely mapped? Can samtools do this? All some other software?
hanleng is offline   Reply With Quote
Old 04-24-2012, 05:00 PM   #2
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Look in the samtools specification... pretty sure there is some type of flag for multiply mapped reads. Then you can convert the file to SAM format and grep -v it.
Heisman is offline   Reply With Quote
Old 04-24-2012, 05:28 PM   #3
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Please define what uniquely mapped means to you.
nilshomer is offline   Reply With Quote
Old 04-24-2012, 05:39 PM   #4
hanleng
Member
 
Location: Houston, TX

Join Date: Mar 2012
Posts: 15
Default

Quote:
Originally Posted by nilshomer View Post
Please define what uniquely mapped means to you.
Reads that can be hit only once in the genome - I know there is some cutoff for this hit, but generally, how to remove those have mutilple hits?
hanleng is offline   Reply With Quote
Old 04-24-2012, 05:40 PM   #5
hanleng
Member
 
Location: Houston, TX

Join Date: Mar 2012
Posts: 15
Default

Quote:
Originally Posted by Heisman View Post
Look in the samtools specification... pretty sure there is some type of flag for multiply mapped reads. Then you can convert the file to SAM format and grep -v it.
That is the reason I asked because I could not find it.
hanleng is offline   Reply With Quote
Old 04-24-2012, 06:04 PM   #6
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

I don't like this definition, since suppose a mapper X tries harder to find a hit than mapper Y, then mapper X is most likely more sensitive and specific, but will have fewer reads that only had one hit, even though a hit for a read may be much more likely than all other hits.

That's why you should look at mapping quality. Those with mapping quality zero are ambiguous: multiple equally best alignments were found. But those with higher quality should have higher confidence.
nilshomer is offline   Reply With Quote
Old 04-25-2012, 12:03 AM   #7
Chipper
Senior Member
 
Location: Sweden

Join Date: Mar 2008
Posts: 324
Default

In other words, use samtools view -q 1 on the .bam to get reads with a mapping quality of at least 1. Depending on your application and aligner you may want to use something like -q 20 to get more reliable hits.
Chipper is offline   Reply With Quote
Old 04-25-2012, 04:38 AM   #8
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

Quote:
Originally Posted by Heisman View Post
Look in the samtools specification... pretty sure there is some type of flag for multiply mapped reads. Then you can convert the file to SAM format and grep -v it.
If the SAM/BAM file has secondary mappings recorded (which not all mapping software will do), then yes, you can filter them out using the FLAG bit values. However I'd recommend using 'samtools view' with the -f and/or -F switches rather than grep.
maubp is offline   Reply With Quote
Old 04-25-2012, 09:23 AM   #9
Simon Anders
Senior Member
 
Location: Heidelberg, Germany

Join Date: Feb 2010
Posts: 994
Default

If your aligner properly sets the MAPQ field, you can filter on this. After all, if there are two equally plausible alignments, the probability for each alignment to be correct is at most 50%, which transforms to a Phred score of 3. Hence, an aligner should never indicate a mapping quality above 3 for multiply matched reads.

Furthermore, many aligners follow the recommendation to use the optional tag "NH" to indicate how many alignments are reported for this read. This helps only, of course, if you instructed the aligner to report multiple alignments.
Simon Anders is offline   Reply With Quote
Old 08-25-2015, 05:04 AM   #10
Annibal
Member
 
Location: Italy

Join Date: Mar 2012
Posts: 10
Default

Quote:
Originally Posted by Simon Anders View Post
If your aligner properly sets the MAPQ field, you can filter on this. After all, if there are two equally plausible alignments, the probability for each alignment to be correct is at most 50%, which transforms to a Phred score of 3. Hence, an aligner should never indicate a mapping quality above 3 for multiply matched reads.
What about paired end sequencing? Uniqueness of a read alignment depends also on its mate position. One mate could map 10 times on the genome but only 1 position is valid considering its mate. How can you filter from a bowtie/bwa generated bam file only uniquely mapped paire end reads?
Don't want considering only "concordant" reads, since i would like to retain paired reads that map in the same region but with discordant orientation (if stranded)
Annibal is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:51 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO