If alignment is done allowing a read to map to multiple locations, is there an easy way to count only the uniquely mapped reads? HTSeq-count "... does not count reads that map to multiple genes.". However, for GAlignments, "... multi-reads (i.e. reads with multiple hits in the reference) won't receive any special treatment ..." meaning that the reads mapping to multiple genomic locations will be counted many times. The reason for allowing a read to map to multiple locations is to make the one of the two resulting BAM files for each sample from STAR useful for RSEM FPKM calculation but also being able to get ordinary counts of uniquely mapped reads for modelling using GLMs.

The STAR manual says that alignments are unique when the NH tag is 1, or when the map quality is 255. One could use these in ScanBamParam(), e.g., ScanBamParam(mapqFilter = 255) or ScanBamParam(tagFilter=list(NH = 1)) or scanBamParam(mapqFilter = 255, tagFilter = list(NH = 1)).

The HTseq documentation seem to suggest that a slightly different heuristic is used, since there is no knowledge of aligner assumed. The heuristic uses a map quality of 10 as a threshold to filter multiple aligners. Also, the documentation implies an 'or' condition on the NH tag, which may or may not be present; ScanBamParam() does not allow this flexibility, but for STAR alignments it seems HTSeq would use the equivalent of ScanBamParam(mapqFilter = 10, tagFilter = list(NH =1)).