How To Tell How Well A Sam File Was Mapped

Question

[Applet]

I'm comparing using STAMPY vs BWA to map some reads to a reference genome.

Now that I've generated two sam files how do I tell how well the mapping did?

Answer 1

[Swbarnes2]

Convert the .sam to .bam with samtools view. .sams are huge. samtools flagstat will take a .bam, and tell you how many reads mapped total.

[Swbarnes2]

I think that's probably true. You could use awk to filter out reads that have mapping quality of 0 from your .sam. I'm not sure that you'd expect a great deal of difference in how many reads have mapping quality of 0 between two aligners, since that's mostly about what your reference sequence is, and how long your reads are.

[David Langenberger]

I fear samtools flagstat counts reads mapping multiple times not unique. Correct me, if I'm wrong.

[David Langenberger]

Well...even a small difference that (might) happen, are differences, right? ;)

[David Langenberger]

Well...even small differences that (might) happen, are differences, right? ;)

[Alby]

Quick Question: Does bam file need to be sorted in order to get the correct flagstat file?

[Sean Davis]

The sorting does not have an effect on flagstat, I do not believe.

[David Langenberger]

I would recommend you to check it. Run flagstat on both files and compare the results.

[Alby]

Just confirmed.BAM file does not need to be sorted for flagstat!

Answer 2

[Sean Davis]

One can generate some pretty useful stats using Picard tools:

http://picard.sourceforge.net/picard-metric-definitions.shtml#AlignmentSummaryMetrics

http://picard.sourceforge.net/picard-metric-definitions.shtml#DuplicationMetrics

If this is a capture experiment, you can also try:

http://picard.sourceforge.net/picard-metric-definitions.shtml#HsMetrics

Answer 3

[David Langenberger]

e.g. you can count the mapped reads and use the percentage of mapped reads as an indicator of how well the mapping tool performed.

Answer 4

[David Langenberger]

In order to assure, you don't have an M somewhere else (read-ID, mitochondrial chromosome, ...), or even reads mapping to multiple positions in the genome (both would basically destroy your statistics), I recommend to do something like this:

perl -F'\t' -ane 'if($F[5]=~/M/){print "$F[0]\n";};' file.sam | sort | uniq | wc

short explanation: if the 6th column contains a 'M' (a match exists), print the first column (the read-ID). then sort and uniq the read-IDs, to assure you count multiple mapped reads only once. and then the word-count as ever.

Answer 5

[Thomasvangurp]

Use grep to find reads that match to the reference. Generally speaking, only lines with matches contain 'M' at position 6 in each line . Use cut -f 6 file.sam| grep M | wc -l

[David Langenberger]

In order to assure, you don't have an M somewhere else (read-ID, mitochondrial chromosome, ...), or even reads mapping to multiple positions in the genome (both would basically destroy your statistics), I recommend to do something like this:

perl -F't' -ane 'if($F[5]=~/M/){print "$F[0]n";};' file.sam | sort | uniq | wc -l

short explanation: if the 6th column contains a 'M' (a match exists), print the first column (the read-ID). then sort and uniq the read-IDs, to assure you count multiple mapped reads only once. and then the word-count as ever.

[David Langenberger]

In order to assure, you don't have an M somewhere else (read-ID, mitochondrial chromosome, ...), or even reads mapping to multiple positions in the genome (both would basically destroy your statistics), I recommend to do something like this:

perl -F'\t' -ane 'if($F[5]=~/M/){print "$F[0]\n";};' file.sam | sort | uniq | wc -l

short explanation: if the 6th column contains a 'M' (a match exists), print the first column (the read-ID). then sort and uniq the read-IDs, to assure you count multiple mapped reads only once. and then the word-count as ever.

[Christof Winter]

That doesn't work if you have an M somewhere else, e.g. in your read name

[Thomasvangurp]

@Christof thanks. Edited the response