I have been using the -q
option of samtools view
to filter out reads whose mapping quality (MAPQ) scores are below a given threshold when mapping reads to a reference assembly with either bwa mem
or minimap2
. However, using samtools idxstats
to count total mapped reads and unmapped reads indicates that these reads with lower MAPQ scores are actually being discarded altogether, since the number of unmapped reads reported (last row) is always zero when the BAM file is created with a -q
argument, even though fewer reads are being mapped to the assembly.
For example:
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b | samtools sort --reference assembly.fasta > nanopore_mapped_no_q_thresh.bam
samtools idxstats nanopore_mapped_no_q_thresh.bam
gives
contig_1 4751250 9897 0
contig_2 37697 528 0
contig_3 9131 2366 0
contig_4 98997 701 0
contig_5 4539 365 0
contig_6 3933 3164 0
* 0 0 142
while
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b -q 2 | samtools sort --reference assembly.fasta > nanopore_mapped_q2.bam
samtools idxstats nanopore_mapped_q2.bam
gives
contig_1 4751250 9716 0
contig_2 37697 528 0
contig_3 9131 1415 0
contig_4 98997 700 0
contig_5 4539 365 0
contig_6 3933 3161 0
* 0 0 0
Is there an alternative way to filter mapped reads by their MAPQ scores so that the filtered reads are still retained in the BAM file?