Tl;dr - can low quality (non-unique) alignments (according to bam_stat.py) be considered multimappers by featureCounts?
I have an RNAseq bam file (subsetted for troubleshooting) that when I run featureCounts, it tells me I have multimapping reads. However, this bam file has been filtered for primary reads only
(-F 256), and my bam_stat output doesn't indicate any multimapped reads:
bam_stat.py -i UWNsub100k_256_noribos.bam Total records: 103011 QC failed: 0 Optical/PCR duplicate: 0 Non primary hits 0 Unmapped reads: 54891 mapq < mapq_cut (non-unique): 11891 mapq >= mapq_cut (unique): 36229 Read-1: 18198 Read-2: 18031 Reads map to '+': 18009 Reads map to '-': 18220 Non-splice reads: 36017 Splice reads: 212 Reads mapped in proper pairs: 32780 Proper-paired reads map to different chrom:0
From the above, it looks like I have 11891 non-unique reads.
When I run featureCounts, unmapped reads and assigned reads appear to agree with the bam_stat output, but featureCounts seems to by considering my non-unique reads as multimappers. Should it not be categorizing them as Unassigned_MappingQuality? And shouldn't multimapping reads be filtered out with
featureCounts -T 4 -p -a ~/path/sacCer3.ncbiRefSeq.gtf -t exon -g gene_id -o ~/path/UWNsub100temp.txt UWNsub100k_256_noribos.bam Assigned 14434 Unassigned_Unmapped 27144 Unassigned_Read_Type 0 Unassigned_Singleton 0 Unassigned_MappingQuality 0 Unassigned_Chimera 0 Unassigned_FragmentLength 0 Unassigned_Duplicate 0 Unassigned_MultiMapping 6034 Unassigned_Secondary 0 Unassigned_NonSplit 0 Unassigned_NoFeatures 3657 Unassigned_Overlapping_Length 0 Unassigned_Ambiguity 281
Any insights are much appreciated - thank you!