1

I want to find and compare the results from STAR and BWA- MEM mapped. I have 150bp paired end reads in sorted.bam files in each case and i want to find in each case uniquely mapped reads, number of reads mapped to multiple loci and number of unmappped reads.

Especially for uniquely mapped reads i tried different solutions from biostars and this Obtaining uniquely mapped reads from BWA mem alignment but i dont know what to trust,

Could you please suggest the best way to find mapping results from bam paired end files ?

  • 1
    Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. – Community Feb 14 '22 at 14:22
  • It would be helpful to describe the genome and a bit of project background regarding what you have already tried. Biostars is a different site to Stackexchange. – M__ Feb 14 '22 at 14:28
  • I have paired end mapped reads from human in sorted bam files and i want to find the unique mapped , the multimapped and the unmapped reads from BWA and from STAR. I already have the .bam files from each mapper. i Tried samtools flagstat and the solution above but i am not sure if this is correct especially for unique reads – user3683485 Feb 14 '22 at 15:08

1 Answers1

2

Bam files contain a MAPQ field that informs on whether an alignment is unique or not. With STAR, a MAPQ value of 255 means the read is mapped uniquely and a score lower than 255 means not unique.

You can use samtools on a bam file to count or output the uniquely or multimapped reads by filtering for the MAPQ as described here: https://groups.google.com/forum/#!topic/rna-star/ipHc9FrgaPU

For example, count uniquely mapped reads: samtools view -c -q 255 myBamFile.bam
Output uniquely mapped reads: samtools view -q 255 myBamFile.bam

There are differences between STAR and BWA MAPQ output though. see here for details: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/