Compensating for multi-mapper count in mpileup?

[Brian.R]

I use mpileup files to call variants for subsections in the genome at a time with a set of custom python/awk scripts. This works well for unique mapping reads, but if I have reads that map to the genome more than once, this information gets lost during mpileup generation.

I have written some scripts to count the 'mappability' of a given read, as my mapper of choice, bowtie, does not report this in the NH or IH tags, but am struggling a bit incorporating this 'correction factor' in to mpileups.

From what I can make out in the manuals samtools might not be able to do so. Am I missing something, or is there perhaps a different pileup-like format I could use to extract 'normalised'-coverage and variants?

[Brian Bushnell]

Depending on your goal, you could simply place multi-mapping reads at a random location, which would accomplish the normalization you are looking for.

[Brian.R]

Well, essentially, if a read maps to 2 loci, and I'm generating a pileup for one of these, I want the bases of this read to only count half towards the total coverage in the pileup.

Does that make sense?

I'm working with miRNAs, so mapping results to multiple identical loci, need to be accounted for in some way (ie dividing by the number of identical loci).

[Brian Bushnell]

Right... that's what random assignment does, on average. A multi-mapping read is only printed as a single line in the sam file, choosing one of the possible mapping locations at random. Many aligners support this as an output option (I think bowtie does).

[Brian.R]

I see what your saying. I was working with those settings beforehand, but discussions with colleagues made me rather uncertain if the bowtie settings reported random positions, or rather just "the first occurrence in the genome".

Probably why I ended up down the path of going for reporting of all mapping positions in the first place.

Thanks.

[Brian Bushnell]

https://www.biostars.org/p/55237/

Looks like there may be some problem in Bowtie's random placement, so I guess you should not rely on it. But since you are postprocessing multimapped reads anyway, YOU could pick one at random and discard the others.

[Brian.R]

Thanks for that link. I'll give it a go with bwa tomorrow and see how that works out.

[Brian Bushnell]

Thanks for that link. I'll give it a go with bwa tomorrow and see how that works out.

Ok. Though I would feel remiss if I did not point you toward BBMap

It supports random assignment for multi-mapping reads with the "ambig=random" flag.

[Brian.R]

Heh, cool, thanks for the suggestion. I'll compare the three of them once I have some time next week or so