how to eliminate reads that map to multiple places in genome

[ammoon]

I am trying to compare expression of histone genes in an rnaseq dataset using igv; many of the reads map to multiple places in the genome,,not surprising given the homology in this protein family but, i need to get accurate expression data . How do i fix the alignments to get rid of these kinds of reads?

[Brian Bushnell]

Typically, ambiguously-mapped reads are either given a single random mapping location, or mapped to all locations via secondary alignments. Expression quantification tools understand this.

You can get rid of multimapping reads by filtering out everything with a MAPQ of 3 or less, but that will cause underrepresentation of the expression of homologous genes.

[ammoon]

thanks for the response

i tried setting the igv threshold for map quality to 50 and i am still getting reads that map to multiple places based on the "blat read sequence" readout table

some reads start in one histone gene and end in another one 90Kb away...?

how do i tell if my reads were given a single random location or mapped to all via secondary alignments?
thank you

[GenoMax]

It depends on what aligner you used and if it followed the consensus of giving a low MAPQ value to the multi-mapped reads. Brian was referring to default behavior of his BBMap aligner. You can either re-map the data (if you are not sure what your original aligner did) or use one of the suggestions in this Biostars thread to filter your BAM files.

[Brian Bushnell]

You can also look at the NH tag:

NH i Number of reported alignments that contains the query in the current record

That's an optional field, but it should tell you how many alignments there are for a read, in the sam file. There's also

H0 i Number of perfect hits

...but that is much less common.

[aprice67]

Check out figure 12 and the area around it in this paper.

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0180904

It might not be such a good idea to remove them.

[Markiyan]

Some commandline examples for mapping quality filtering.

Mapping score distribution depends on mapper program/version used.

First check the actual MAPQ score histogram for your bam file (first 1M of reads):

samtools view $1 | head - -n 1000000 | cut -f5 | sort | uniq -c

Than filter with desired threshold using samtools and -q parameter (assuming the input is test.bam):

#!/bin/sh
MAP_QUAL=30
STRAIN=test
samtools view -q $MAP_QUAL -b ${STRAIN}.bam -o ${STRAIN}.Q$MAP_QUAL.bam -U ${STRAIN}.Q${MAP_QUAL}U.bam


The passing MAPQ filter reads will end up in test.Q30.bam

And failing ones would end up in test.Q30U.bam file.

[aprice67]

Be very careful when dealing with MAPQ scores. It seems every aligner has a different way to determine them and they aren't standardized at all. See here:

http://www.acgt.me/blog/2015/3/17/more-madness-with-mapq-scores-aka-why-bioinformaticians-hate-poor-and-incomplete-software-documentation