2
$\begingroup$

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?

$\endgroup$

1 Answer 1

2
$\begingroup$

Looking at the samtools view docs, I believe that you need to also set the -p flag if you want to retain "filtered" reads in the file and simply have them be unmapped:

The -p, option sets the UNMAP flag on filtered alignments then writes them to the output file.

$\endgroup$
1
  • 2
    $\begingroup$ Yes, this does the trick. This flag was not recognized by the older version of samtools I was using (1.11) but after upgrading to 1.15 it works. Thank you Maximilian! $\endgroup$ Jul 24, 2022 at 3:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.