SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
TopHat/Bowtie - number of reads aligned mgibson Bioinformatics 7 10-22-2011 09:04 PM
How many uniquly aligned sequences can we expect? chubbchubb RNA Sequencing 0 02-11-2011 09:14 AM
Extending aligned sequences, plus/minus strand biznatch Bioinformatics 3 01-21-2011 11:31 AM
Determining the number of bases and percent coverage in an aligned sequence kz26 Bioinformatics 0 06-28-2010 09:22 PM
the number of aligned reads by maq and bwa totalnew Bioinformatics 3 12-14-2009 12:18 PM

Reply
 
Thread Tools
Old 06-10-2010, 06:19 AM   #1
Cbon
Junior Member
 
Location: Paris, France

Join Date: Jun 2010
Posts: 2
Default BWA - number of sequences aligned

Hello everyone,

I am quite a newbie in bioinformatics, and I have some difficulties using bwa software.

I successfully aligned my dataset to my reference genome. Now, I want to see the impact of the settings on my results, but I don't know how to obtain the number of reads aligned in the .sam file.

Is there a useful trick to obtain the number of reads aligned ?

I have found no clue on the bwa manual, neither in the samtool package and I begin to feel quite helpless ...

Thank you for your answers !
Cbon is offline   Reply With Quote
Old 06-10-2010, 06:39 AM   #2
dawe
Senior Member
 
Location: 4530'25.22"N / 915'53.00"E

Join Date: Apr 2009
Posts: 258
Default

Quote:
Originally Posted by Cbon View Post
Hello everyone,

Is there a useful trick to obtain the number of reads aligned ?
Try

Code:
samtools flagstat file.bam
Although I don't know if it works on sam files (in case you should convert in BAM, which I believe should be the default...).

d
dawe is offline   Reply With Quote
Old 06-10-2010, 06:44 AM   #3
racetrout
Junior Member
 
Location: North Pole

Join Date: Jun 2010
Posts: 2
Default

If you are using UNIX and GNU awk,
Code:
gawk '!/^[@#]/ && (!and($2,4){print}' bwa.sam| wc -l
should give you the number of reads aligned - providing BWA is setting the flag for unaligned reads. If you have many reads and BWA writes out only the aligned reads,
Code:
wc -l bwa.sam
is a very good simple approximate of aligned reads (it counts the header of course, but the ~5 lines of the header compared to the millions of reads are negligible).
racetrout is offline   Reply With Quote
Old 06-10-2010, 07:04 AM   #4
Cbon
Junior Member
 
Location: Paris, France

Join Date: Jun 2010
Posts: 2
Default

Thank you both for your answers !

My files are already converted in .bam, so I hope I can use the flagstat command. If it doesn't work, I will try the code.

thank you again !
Cbon is offline   Reply With Quote
Old 05-13-2011, 08:28 AM   #5
saml
Junior Member
 
Location: Uppsala, Sweden

Join Date: Oct 2010
Posts: 4
Default

Quote:
Originally Posted by racetrout View Post
If you are using UNIX and GNU awk,
Code:
gawk '!/^[@#]/ && (!and($2,4){print}' bwa.sam| wc -l
should give you the number of reads aligned - providing BWA is setting the flag for unaligned reads. If you have many reads and BWA writes out only the aligned reads,
Code:
wc -l bwa.sam
is a very good simple approximate of aligned reads (it counts the header of course, but the ~5 lines of the header compared to the millions of reads are negligible).
Yes, but removing the lines is easy too. Skipping 2 header lines (think I had that), and starting on the 3rd line, can be done with the tail command and adding a '+' before the number of lines, to denote which line to start from:

Code:
tail -n +3 bwa.sam | gawk '!/^[@#]/ && (!and($2,4){print}' | wc -l
... or for the case with only aligned reads coming from bwa:
Code:
tail -n +3 bwa.sam | wc -l
saml is offline   Reply With Quote
Old 05-13-2011, 09:58 AM   #6
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

bwa by default includes all reads, even unaligned ones. You can tell that they are unmapped, from the binary tag.

So the flagstat command will distinguish between the mapped an unmapped reads.

You can also use samtools view -bF0x4 all.bam > mapped .bam to get a bam which has all the unmapped reads filtered out.
swbarnes2 is offline   Reply With Quote
Old 05-13-2011, 10:24 AM   #7
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Use "samtools view -c -S <in.sam>" and add the "-f/-F" options as desired.
nilshomer is offline   Reply With Quote
Old 05-13-2012, 04:45 PM   #8
solonas13
Junior Member
 
Location: Germany

Join Date: May 2012
Posts: 1
Default

I personally find it---at the best case---unacceptable the fact that a short-read aligner does not have an extra line of code to report the total number of aligned reads... Especially when it outputs (on the standard error) things like:

`[infer_isize] (25, 50, 75) percentile: (137, 179, 227)'

Anyhow, isn't there a closing parenthesis missing from this one?

Quote:
Originally Posted by saml View Post
Code:
tail -n +3 bwa.sam | gawk '!/^[@#]/ && (!and($2,4){print}' | wc -l
Thanks

Last edited by solonas13; 05-13-2012 at 04:47 PM.
solonas13 is offline   Reply With Quote
Reply

Tags
bwa, sam, samtools

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:38 PM.


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