Using samtools to read a bam file

Question

[Caitlin]

Hi all.

I am trying to extract some information (please see below) from a bam file using samtools. Given my lack of experience, I am struggling.

1.) Total reads aligned on the forward strand 2.) Total reads aligned on the reverse strand 3.) Number of read pairs where the pairs are pointing towards each other 4.) Number of read pairs where the pairs pointing in the same direction (this is the opposite of proper pairs) 5.) Number of primary alignments (not SECONDARY) that are on the forward strand

Any help would be greatly appreciated.

Thank you.

[Chris Dean]

You can use the samtools stats command to extract some of this information (version samtools-1.13). More information about this command can be found on the manual page. Hope this helps!

[Caitlin]

Thank you Chris. Unfortunately, the only version of samtools I have access to is version 1.9 I ran samtools stats on the bam file, but I could not locate any of the information I needed.

[Chris Dean]

Interesting, I just ran samtools (v1.9) on a BAM file and I do see some of this information. For example, are the 1st fragments, last fragments, and inward oriented pairs counts not what you are looking for? These outputs can be found in the SN section of the output.

[Caitlin]

Thanks Chris. I'll look again.

1st fragments = Total reads aligned on the forward strand last fragments = Total reads aligned on the reverse strand inward oriented pairs = Number of read pairs where the pairs are pointing towards each other

??

I was able to upgrade samtools to 1.12.0 which should, hopefully, be of some use.

[Chris Dean]

Yes exactly.

[Caitlin]

Is there a field for the last two questions? Number of read pairs where the pairs pointing in the same direction (this is the opposite of proper pairs) and Number of primary alignments (not SECONDARY) that are on the forward strand?

[Caitlin]

Hi Chris.

I thought I would include the link to the file I am using.

http://data.biostarhandbook.com/bam/demo.bam

Does this help? If you could run the command and compare answers with me that would be ideal.

Thank you Chris.

Answer 1

[swbarnes2]

Why don't the binary flags give you exactly this information?

[jkbonfield]

The flag field only distinguishes "properly paired" or not. The definition of that is up to the aligner.

It's not even that obvious either. Some old mate pair libraries has large inserts with reads pointing away from each other, where the enrichment for circularised templates worked, or short inserts with the pairs pointing towards one another (where the enrichment failed). Both can be mapped in a consistent manner. I doubt any aligners coped in setting "properly paired" for both portions though!

[swbarnes2]

The direction fields should work, even if for some reason the properly paired designation is wrong. And why should it be wrong for huge numbers of reads? OP hasn't indicated that he is doing anything exotic.

[Caitlin]

Hi.

If you're referring to me, I'm a female.