Multi-mapped reads when mapping to reference genome

[michael-kotliar]

I believe that "$PARAM_BAM_PREFIX"_total.bam should be either filtered by quality 255 like you do here or by flag 0x100 (not primary alignment) similar to the way you do here. Otherwise your BAM file includes multi-mapped reads, that influences on scaling coefficient. Alternatively, you can use --outFilterMultimapNmax 1 for STAR aligner and don't use any filters after it. By default, STAR allows up to 10 multiple alignments.

[julienrichardalbert]

Filtering the reference-aligned bam file with the same criteria for filtering pseudogenome-aligned bams was considered.. in fact it is still an active topic in our group. One pragmatic reason it is coded as such is because we are interested in multi-mapped reads (especially to repetitive elements).

The main reason why we keep multi-mapped reads is because "mappability" is different between the reference and pseudo-genomes. This is important at genes part of large families - when a SNV or INDEL is introduced, this gene will become "more mappable" in the pseudogenome relative to the reference. Then, when creating the genomic tracks, if uniquely aligned reads are removed from the reference bam, allele-specific reads outnumber "total" reads at these repetitive genes. This causes really bad visualization of allelic effects, which is something we want to avoid.

A word of caution in specifying the number of multimap alignments to report: when using Tophat2, setting -g/--max-multihits resulted in unique reads alignments at repetitive loci. Instead of discarding all alignments with equal mapping quality, the aligner marks one of the alignments as unique (high MAPQ). Not sure if this applies to STAR as well, but I am now wary of these parameters.

Happy to discuss more on this issue