I would like to get all the alignment positions of Illumina reads on several closely related genomes (if the aligments are above a certain aligment threshold).
So I use the the "-a" flag to output all the aligments. The problem is, that those reads get flaged as secondary.
How do I distinguish between secandary hits which are split hits (So I have to use "-M" flag as well) and the aligments which are marked secondary because a read hits multiple positions on the genomes with same/similar qual.
A small sidequestion: If a read has multiple hits on different genomes with the same mapping quality, the hit which gets reported is randomly determined, right?
Not answering your question directly but you could use
bbmap.sh
with multiple genome references and use the optionambig=all
(in addition to any other flags you need) to allow reads to multi-map to all genomes.