Is there a way to do read filtering (MAPQ> certain value) on a BAM instead of SAM file?

Question

[msimmer92]

I'm working on Bash scripts for a ChIPseq pipeline for my lab. Even though the ENCODE guideline suggests to remove duplicates, some people here want to not remove duplicates but filter the reads with certain MAPQ values. For this purpose, I am working on a script that does that.

On the Internet and some other people's scripts, all the examples I have seen so far are filtering SAM files (for example, with awk or grep, knowing that the MAPQ value is on the 5th column of the SAM file, it's not a challenge to extract this; let's say it becomes a simple file-management and text-editing problem). Nevertheless, in the pipeline I'm working on, the inputs come in the form of sorted BAMS (because there's another script in the lab that does the mapping, sorting and conversion to BAM).

So I was wondering, is there a way of doing this filtering of MAPQ=certain values from the sorted BAMs I got from the people, without having to ask them for the SAM files? Thank you!

[msimmer92]

Thank you, now it's all clear :) . One more question: if instead of wanting to filtering a BAM for mapping quality I wanted to filter out reads with q<20 in a fastq file (so fastq filtering, not bam), which would be the most recommendable way to do it? Thank you again

[h.mon]

Fastq files contain quality score per base, maybe you want reads with average quality < 20?

[msimmer92]

Considering both possibilities - per base filtering and per average read score - what is your recommendation? The concept is that I would like to have a fastq file of higher quality, but I'm not sure which approach is best to take.

Answer 1

[h.mon]

-q INT only include reads with mapping quality >= INT [0]

    samtools view -h -q 20 file.bam

[msimmer92]

This is great, thank you ! And if I wanted to include two conditions? (such as include reads that MAPQ==1 and MAPQ>20 )?

[Bastien Hervé]

Why would you want to keep MAPQ=1 ? This is pretty bad mapping quality. What are your trying to do ?

[msimmer92]

A labmate has a script that does that in order to also keep uniquely-mapped reads (which if you have a SAM file made by Bowtie2, it's MAPQ==1), because it seems it's good to have them in ChIPseq data analysis. I'm looking at his code to see if his has certain functions (like that, for example) that mine doesn't and might me good to implement (so, I am open to all your opinions if you think that is not correct/necessary). We always use Bowtie2 for ChIPseq mapping, that's why we use those numbers in particular (uniquely mapped reads ==1 and reads with quality >4). I clarify because depending on the mapper, these numbers might change.

[Bastien Hervé]

MAPQ is a mapping quality, not the number of multi mapped location for a read.

If you want to keep uniquely-mapped reads + MAPQ > 20

samtools view -h -F 256 -q 20 file.bam

[msimmer92]

Thank you, I will show this to the people in my lab. I agree with you. There is a mistake in his code, then. I read also the Bowtie2 manual in detail and also states the higher the MAPQ the higher the "uniqueness" of the read, so it doesn't make sense to keep the MAPQ==1 reads.

[Bastien Hervé]

Actually, I never filter mapped reads on mapping quality. If a read is badly mapped on the reference, Bowtie2 will not align it. I say that, the mapped reads fit good enought to be mapped so I keep them all.

Also, take a look at this

[msimmer92]

To be clear, you use "-F 256", that is the flag of the secondary-alingment reads, to filter those ones in particular, right? (I'm making sure I got that right). The logic would be "from the secondary alignment reads, filter out those ones with q smaller than 20"

[Bastien Hervé]

No, with this command you filter out not primary alignment and you also filter out primary alignment with less than 20 MAPQ. In other words you keep only primary alignment with MAPQ >= 20

[msimmer92]

Thank you, it's clear now. One clarification: is this depending on the sorting of the bam? (If I do this filtering without having the bam sorted, would I have a problem? ) Thank you

[Bastien Hervé]

I don't think sorted bam is necessary to filter out but it's always a good point to sort it