flagstat

[jwhite]

Hi Folks,

Our standard process for preparing BAM files involves the following:

samtools sort
samtools rmdup
samtools sort
samtools index

20682963562 Mar 19 20:13 EAP120_R2.fastq.bam
12717567745 Mar 19 21:50 EAP120_R2.fastq.srt.bam
4529263043 Mar 19 22:11 EAP120_R2.fastq.nodup.bam
4529263062 Mar 19 22:38 EAP120_R2.fastq.nodup.srt.bam

So there will be 4 BAM files from this process. To evaluate them, I used samtools flagstat.
What I found with one of our samples is that the number of mapped reads drops from 17.8 billion to 5.8 billion. (The flagstat data are below.)

Why would that occur? Does that indicate poor quality data?
It would seem to me to indicate that there are a LOT of duplicates. Is that correct? Should I normally see a 75 % reduction in the number of mapped reads after duplicate removal?

So, I guess what I'm asking is how should I interpret flagstat results?

Thanks for any light you can shed on this.

Joe White

-----------------------
# initial BAM file
samtools flagstat EAP120_R2.fastq.bam
177833072 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
172622167 + 0 mapped (97.07%:-nan%)
177833072 + 0 paired in sequencing
88916536 + 0 read1
88916536 + 0 read2
169943152 + 0 properly paired (95.56%:-nan%)
171556513 + 0 with itself and mate mapped
1065654 + 0 singletons (0.60%:-nan%)
1376256 + 0 with mate mapped to a different chr
1254219 + 0 with mate mapped to a different chr (mapQ>=5)

# after the first sort
samtools flagstat EAP120_R2.fastq.srt.bam
177833072 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
172622167 + 0 mapped (97.07%:-nan%)
177833072 + 0 paired in sequencing
88916536 + 0 read1
88916536 + 0 read2
169943152 + 0 properly paired (95.56%:-nan%)
171556513 + 0 with itself and mate mapped
1065654 + 0 singletons (0.60%:-nan%)
1376256 + 0 with mate mapped to a different chr
1254219 + 0 with mate mapped to a different chr (mapQ>=5)

# after duplicate removal
samtools flagstat EAP120_R2.fastq.nodup.bam &
57792684 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
52581779 + 0 mapped (90.98%:-nan%)
57792684 + 0 paired in sequencing
28898595 + 0 read1
28894089 + 0 read2
50015042 + 0 properly paired (86.54%:-nan%)
51516125 + 0 with itself and mate mapped
1065654 + 0 singletons (1.84%:-nan%)
1376256 + 0 with mate mapped to a different chr
1254219 + 0 with mate mapped to a different chr (mapQ>=5)

# after the second sort
samtools flagstat EAP120_R2.fastq.nodup.srt.bam &
57792684 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates
52581779 + 0 mapped (90.98%:-nan%)
57792684 + 0 paired in sequencing
28898595 + 0 read1
28894089 + 0 read2
50015042 + 0 properly paired (86.54%:-nan%)
51516125 + 0 with itself and mate mapped
1065654 + 0 singletons (1.84%:-nan%)
1376256 + 0 with mate mapped to a different chr
1254219 + 0 with mate mapped to a different chr (mapQ>=5)

[swbarnes2]

You don't need to sort after the rmdup.

But yes, sometimes a library really is that repetitive. Maybe you have way more reads than your small genome needs. Maybe your library prep got screwed up.

Rule #1 Don't disbelieve your data just because you don't like what its telling you.

If you like, use Picard MarkDuplicates instead, because it will only mark, and not remove the duplicates, and then you can eyeball and see the repetition. Or look at the first sort file in IGV, or something like that, and see if you can see the duplicates.

[mmmm]

using flagstat- what does this mean?
"76 reads have mapped to a different chromosome than their pair" (BWA was used for reads mapping to the reference bacterial genome)

[dpryan]

The origin of that is comparison of the RNAME and RNEXT columns. If they're not the same, then this is incremented. I would have to presume that the bacterial genome you aligned against contained not just the main genome but also some plasmids (or perhaps it was all just a bunch of contigs).

[mmmm]

yes- there is a plasmid but it showed only 1 SNP

I have got these stats for mapping to the chromosome and plasmid- does "76 + 0 with mate mapped to a different chr" refere to repitative sequences that are shared between plasmid and chromosome

3753203 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
3660519 + 0 mapped (97.53%:-nan%)
3753203 + 0 paired in sequencing
1876530 + 0 read1
1876673 + 0 read2
3644736 + 0 properly paired (97.11%:-nan%)
3653854 + 0 with itself and mate mapped
6665 + 0 singletons (0.18%:-nan%)
76 + 0 with mate mapped to a different chr
74 + 0 with mate mapped to a different chr (mapQ>=5)

[dpryan]

Given that 74/76 have MAPQ>=5 it would seem that these are unlikely to be that repetitive. However, they might be technical artefacts, since there are so few of them.

[mmmm]

thanks- if I want to determine the (0.18%) of unmapped reads- what is the simplest way for this?- should I use velvet for denovo assembly (kmer 31) then annotate the output?

[dpryan]

~2.5% and yeah that'd be one approach. You might have a look at the reads first to see if they just didn't align due to low quality or something like that.

[mmmm]

thank you- what do you mean by looking at the reads (do you mean viewing the bam file using IGV)

[dpryan]

They're unmapped, so you would just "zcat unmapped.fq.gz | less" or even just run the unmapped reads through fastqc. The general idea is to first see if they're even worth assembling (if they're all repetitive regions or low quality then there's no reason to bother).