Multimapped and unique mapped for BAM file of RRBS data

Question

[sata72]

Hi all,

I have a BAM files from RRBS data, and, I want to have number of all reads, number of multimapped and number of unique mapped. I think it can be done with samtools, and, with this we can have all reads.

Thanks in advance.

 samtools view -c sample.bam

[GenoMax]

To get you started. Multiple other threads exist.

Multi-mapped reads: How to filter/view reads mapped to multiple locations using samtools or pysam
Unique reads: Help with Samtools flags to get only primary, unique, non duplicate aligned reads

[ATpoint]

The exact definition of what is "unique" depends on the aligner. Typically, one sets a MAPQ threshold, for example 20, and everything below is ignored during analysis because mapping quality is too low. You could simply count reads with MAPQ > 20 (samtools view -q 20) and compare to all reads.

Answer 1

[young_bioinformatician]

This command is not enough to do all you want. As far as I know 0x100 is that reads that are marked as multimapped. So the following commands will give you everything you need: All reads :

samtools view -c input.bam

Multimapped:

samtools view -F 0x4 -f 0x100 input.bam | wc -l

unique map reads;

bedtools bamtobed -i input.bam > output.bed
sort -k 4,4 output.bed | uniq -f 3 -u | wc -l

Then you can merge all the information into one.