Hello,
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 -F 256?
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!
