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!
1 answer
Why would you want to keep MAPQ=1 ? This is pretty bad mapping quality. What are your trying to do ?
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.
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
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.
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
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"
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
I don't think sorted bam is necessary to filter out but it's always a good point to sort it
Log in to answer this question.
-
umitools dedup and quality filtering
written by song •
-
BAI files after removing duplicates and multimappers
written by Maka •
-
bowtie2 MAPQ scores more than 42
written by rakusmaximus •
-
In the NGS pipeline, why read are sorted before marking duplicates?
written by ManuelDB
-
Can BWA-MEM MAPQ Score = 70?
written by Saha.Neiladri •
-
Snakemake SyntaxError: No rule keywords allowed after run/shell/script/wrapper/cwl in rule
written by msimmer92
-
Pindel and GATK AD Concordance
written by jsneaththompson
-
Filtering by FPKM, opinions and thoughts
written by Biogeek
-
Is there a correct order in which to filter BAM files?
written by James Ashmore
-
Samtools And Flagstat
written by nbvasani
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
Fastq files contain quality score per base, maybe you want reads with average quality < 20?
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.