# How to Filter BAM Files by Mapping Quality, Flags, and Regions By the end of this walkthrough you will be able to filter a BAM file by mapping quality, SAM flags, and genomic regions, then check that the output is correct before you use it in variant calling, read counting, visualization, or another downstream step. You will know how to combine the three filters in one `samtools view` command, calculate flag masks without guessing, and catch the silent errors that produce plausible but wrong BAM files. You need an input BAM, a coordinate-sorted copy of that BAM, its index, and a BED file or region string for the target intervals. You also need `samtools` installed and available on your path. `bedtools` is optional for independent region checks. Keep the original BAM unchanged and write every filtered result to a new file. BAM files are binary, compressed representations of SAM alignment records. They are central to next-generation sequencing pipelines, and tools that process them must handle large volumes of alignment data efficiently [1]. That practical detail matters because filtering is often the first step that reduces a large BAM to a smaller, task-specific file. ## Filter BAM Files by Mapping Quality, Flags, and Regions: What Each Filter Means The three filters answer different questions. Mapping quality asks how confident the aligner is about the read's placement. Flags ask what kind of alignment record you are looking at, such as unmapped, duplicate, secondary, or proper pair. Regions ask where the alignment sits on the reference. You can apply all three at once, but each one has its own coordinate rules and its own failure modes. ### Bioinformatics file formats around a BAM filter A BAM is the binary form of a SAM file. SAM, Sequence Alignment/Map, is a text format with a header section and one tab-delimited line per alignment record. BAM stores the same information in a compressed binary layout that is faster for tools to read and smaller on disk. CRAM is another aligned-read format that uses the reference sequence for extra compression. FASTQ and FASTA sit earlier in the workflow. FASTQ stores raw sequence reads plus a per-base quality score for each base. FASTA stores sequence letters without per-base quality scores. Neither FASTQ nor FASTA has mapping quality because mapping quality exists only after reads have been aligned to a reference. A BAM record carries the read sequence, the alignment position, the mapping quality, the SAM flag, the CIGAR string, and other alignment fields. This distinction prevents a common category error. You can trim or quality-filter a FASTQ by per-base quality. You cannot filter a FASTQ by MAPQ. You filter a BAM by MAPQ because the BAM is the file that contains alignment placements. The same BAM can be converted back to FASTQ for some workflows, but the MAPQ value is lost when you go back to raw reads. BAM processing tools often stream records or split work across threads because the files can be large [1]. The commands below use streaming filters, so they do not load the whole BAM into memory. ### Mapping quality is a placement confidence score Mapping quality, or MAPQ, is a per-alignment score that summarizes how confident the aligner is that the read is placed at the reported reference position. In the SAM specification, MAPQ is an integer from 0 to 255. Higher values indicate greater confidence. A MAPQ of 0 usually means the read maps equally well to multiple places, or that the aligner could not assign a confident placement. A MAPQ of 255 is reserved for cases where mapping quality is unavailable. MAPQ is not base quality. Base quality is the per-base probability that a base call is wrong. Mapping quality is the per-read probability that the alignment position is wrong. A read can have high base quality and low mapping quality, especially in repeats, paralogous regions, or multi-copy genes. MAPQ scales are not identical across aligners. BWA-MEM, Bowtie2, minimap2, and other aligners assign MAPQ values with different models and different maximum values. A cutoff of 30 is a common starting point in short-read DNA work, but it is not a universal rule. Check the aligner documentation and inspect the MAPQ distribution in your own BAM before choosing a threshold. ### SAM flags are bitwise labels A SAM flag is an integer that encodes several yes or no properties for an alignment record. Each property is a bit. The flag value is the sum of the bits that are set. For example, a value of 4 means the read is unmapped. A value of 1024 means the read is a duplicate. A value of 1028 means the read is both unmapped and a duplicate because 4 + 1024 = 1028. The SAM format specification defines the standard bit meanings. The table below lists the flags you will use most often when filtering BAM files. | Bit | Hex | Decimal | Meaning | |--|--:|--:|--| | 0 | 0x1 | 1 | Read is paired | | 1 | 0x2 | 2 | Proper pair | | 2 | 0x4 | 4 | Read is unmapped | | 3 | 0x8 | 8 | Mate is unmapped | | 4 | 0x10 | 16 | Read is on the reverse strand | | 5 | 0x20 | 32 | Mate is on the reverse strand | | 6 | 0x40 | 64 | First read in a pair | | 7 | 0x80 | 128 | Second read in a pair | | 8 | 0x100 | 256 | Secondary alignment | | 9 | 0x200 | 512 | Read fails vendor quality checks | | 10 | 0x400 | 1024 | Duplicate read | | 11 | 0x800 | 2048 | Supplementary alignment | In `samtools view`, `-f` requires all bits in the supplied mask to be present. `-F` excludes any record that has any bit in the supplied mask set. This difference is the source of many wrong filters. If you want to keep reads that are paired, `-f 0x1` requires the paired bit. If you want to remove duplicates, `-F 0x400` removes records with the duplicate bit. If you want to remove unmapped, secondary, QC-failed, duplicate, and supplementary records, you build a combined exclusion mask. ### Regions use two different coordinate systems A region filter restricts output to alignments that overlap one or more genomic intervals. You can give `samtools view` a region string such as `chr1:1001-2000`, or you can give it a BED file with `-L regions.bed`. The two forms use different coordinate conventions. A region string is 1-based and inclusive. `chr1:1001-2000` includes positions 1001 through 2000. A BED interval is 0-based and half-open. A BED line with `chr1 1000 2000` covers positions 1001 through 2000 in 1-based coordinates. The same interval can be written either way. Mixing the conventions shifts your target by one base, which is easy to miss because the command still runs. Region filtering with `-L` includes reads that overlap the interval, not only reads that are fully contained inside it. A read can start before the interval and still be returned if part of its alignment overlaps the target. For paired-end data, one mate can overlap the region while the other mate sits outside it. That behavior is normal for alignment-level filtering. ## Step-by-Step BAM Filtering Workflow The workflow below assumes `input.bam` is your starting file. The examples use `input.sorted.bam` after sorting and indexing. Replace names with your own file names. Run the commands in a shell where `samtools` is available. ### Step 1: Check that the BAM is readable and indexed Start with file integrity and basic counts. `samtools quickcheck -v` checks whether the file has a valid header and end-of-file marker. It prints nothing when the file passes. `samtools view -H` prints the header. `samtools idxstats` prints per-reference counts and requires an index. ```bash samtools version samtools quickcheck -v input.bam samtools view -H input.bam | head samtools idxstats input.bam | head samtools view -c input.bam ``` Representative output from the header and index checks looks like this. Your reference names and counts will differ. ```text @HD VN:1.6 SO:coordinate @SQ SN:chr1 LN:248956422 @SQ SN:chr2 LN:242193529 chr1 248956422 482193 1204 chr2 242193529 391002 983 1000000 ``` The final number is the alignment record count, not the number of read pairs. A paired-end library can have two alignment records for each DNA fragment. Keep that distinction in mind when you compare before and after counts. If `samtools quickcheck` prints an error, stop and repair or re-download the BAM. If `samtools idxstats` reports that it cannot retrieve an index, create one after sorting. ### Step 2: Sort and index if needed Region filtering and many downstream tools expect a coordinate-sorted BAM with an index. If your BAM is already coordinate-sorted and has a fresh index, skip this step. If it is not sorted, sort it to a new file and index the sorted output. ```bash samtools sort -@ 4 -o input.sorted.bam input.bam samtools index -@ 4 input.sorted.bam samtools quickcheck -v input.sorted.bam ``` No output from `samtools quickcheck -v` means the sorted BAM passed the quick check. `samtools index` creates an index file next to the BAM, usually with `.bai` appended. The `-@ 4` option uses four threads for compression or decompression. Increase it if your machine has capacity, but do not set it higher than the number of available cores. If you try to index an unsorted BAM, you may see a message like this. ```text [E::hts_idx_push] Unsorted positions on sequence #1: 12345 followed by 12000 ``` The fix is to sort first, then index the sorted file. Do not index the original unsorted BAM and then sort it. ### Step 3: Profile MAPQ before choosing a cutoff Do not choose a MAPQ cutoff from habit alone. Count how many records pass at several thresholds, then relate that count to your experiment. The loop below reports counts at MAPQ 0, 10, 20, 30, and 60. ```bash for q in 0 10 20 30 60 do printf 'MAPQ >= %s: ' "$q" samtools view -c -q "$q" input.sorted.bam done ``` Representative output looks like this. ```text MAPQ >= 0: 1000000 MAPQ >= 10: 902144 MAPQ >= 20: 861003 MAPQ >= 30: 742318 MAPQ >= 60: 501776 ``` The exact numbers do not matter for the workflow. The shape does. A steep drop between 20 and 30 can mean many reads sit in repetitive regions. A large fraction at MAPQ 0 can mean a multi-mapping library, a contaminated reference, or an aligner with a different MAPQ scale. If the library is RNA-seq, a low MAPQ fraction may reflect genuine multi-mapping transcripts. If the library is short-read whole-genome sequencing, a large MAPQ 0 fraction deserves investigation. A MAPQ cutoff of 20 or 30 is common in short-read DNA workflows. Use it as a starting point, not as a law. Some analyses need a permissive cutoff to keep reads in repeats. Others need a strict cutoff to reduce false positives. Record the cutoff and the reason in your analysis notes. ### Step 4: Calculate the flag mask Suppose you want to exclude unmapped reads, secondary alignments, QC-failed reads, duplicate reads, and supplementary alignments. The relevant bits are 0x4, 0x100, 0x200, 0x400, and 0x800. Add them once, then use the resulting mask in every command. ```bash printf '%d\n' $((0x4 + 0x100 + 0x200 + 0x400 + 0x800)) printf '0x%X\n' $((0x4 + 0x100 + 0x200 + 0x400 + 0x800)) ``` Expected output: ```text 3844 0xF04 ``` So the exclusion option is `-F 0xF04`. The decimal form `-F 3844` is equivalent, but hexadecimal is easier to audit because each digit maps to a known flag block. If duplicates should stay in the BAM, remove 0x400 from the sum. ```bash printf '0x%X\n' $((0x4 + 0x100 + 0x200 + 0x800)) ``` Expected output: ```text 0x904 ``` That mask excludes unmapped, secondary, QC-failed, and supplementary records but keeps duplicates. Use it when duplicate marking has not been done yet or when duplicate reads are part of the analysis. Be careful with `-f 0x2`. That option requires the proper-pair bit. It is not a general paired-end filter. A read can be paired and correctly flagged as paired without being part of a proper pair. Requiring `0x2` can remove discordant pairs, chimeric reads, and reads in rearranged regions. Use it only when the analysis genuinely needs proper pairs. ### Step 5: Build and verify a BED file A BED file for `samtools view -L` needs chromosome names that match the BAM header exactly. If the BAM uses `chr1` and the BED uses `1`, the command can run and return zero reads. That is a silent failure. Create a small BED file with tabs, not spaces. The `printf` command below writes two intervals. `cat -A` shows tabs as `^I` and line ends as `$`. ```bash printf 'chr1\t1000\t2000\nchr2\t50000\t51000\n' > regions.bed cat -A regions.bed ``` Expected output: ```text chr1^I1000^I2000$ chr2^I50000^I51000$ ``` The first interval covers 1-based positions 1001 through 2000. The second covers 1-based positions 50001 through 51000. If you prefer a region string for a single interval, `chr1:1001-2000` selects the same first interval on an indexed BAM. The BED form is easier to reuse for many intervals. Check that the chromosome names exist in the BAM header before you trust a zero count. ```bash samtools view -H input.sorted.bam | grep '^@SQ' ``` If the BED uses `chr1` and the header uses `1`, fix the BED. Do not assume the region filter failed because there are no reads. Check the names first. ### Step 6: Run a dry count with all three filters Run the combined filter with `-c` before you write a BAM. This gives you a count to compare with the final output and catches syntax errors early. The command below keeps records with MAPQ 30 or higher, excludes the five flag classes, and restricts to the BED intervals. ```bash samtools view -c -q 30 -F 0xF04 -L regions.bed input.sorted.bam ``` Representative output: ```text 742318 ``` The count is the number of alignment records that pass all three filters. It is not the number of reads, not the number of read pairs, and not the number of fragments. A paired-end fragment can contribute two passing records if both mates pass. One mate can pass while the other fails. If the count is zero, check the BED chromosome names, the MAPQ threshold, and the flag mask before changing the command. A zero count from a valid BAM usually means a naming mismatch, an interval in the wrong coordinate system, or a flag filter that removes everything. This dry count also avoids writing a large intermediate file just to discover a typo. Filtering a large BAM in one streaming pass keeps memory use predictable [1]. ### Step 7: Write the filtered BAM Write the filtered BAM with `-b` and `-o`. The output retains the BAM header. The `-@ 4` option uses four threads. The input must be indexed for efficient region access. ```bash samtools view -@ 4 -b -q 30 -F 0xF04 -L regions.bed -o filtered.bam input.sorted.bam samtools index -@ 4 filtered.bam samtools quickcheck -v filtered.bam ``` No output from `samtools quickcheck -v` means the new BAM passed the quick check. `samtools index` creates `filtered.bam.bai`. The order matters. Index the filtered BAM after it is written, not before. If you forget the index, a later region query may print a message like this. ```text [E::idx_find_and_load] Could not retrieve index file for 'filtered.bam' ``` The fix is `samtools index filtered.bam`. If the BAM is not coordinate-sorted, sort it first. ### Step 8: Validate the output Start with counts. The filtered BAM count should match the dry count from Step 6. The region overlap count should also match if every retained record overlaps a BED interval. ```bash samtools view -c filtered.bam samtools view -c -L regions.bed filtered.bam samtools view -c -q 30 -F 0xF04 -L regions.bed input.sorted.bam ``` Representative output: ```text 742318 742318 742318 ``` The three numbers should agree. If the first and second differ, some records in the filtered BAM do not overlap the BED intervals. That can happen if you wrote the BAM with different options than the count command. If the first and third differ, the write step did not use the same filters as the count step. Next, check the flags you intended to remove. These commands count records that still carry the excluded bits. Each should return zero. ```bash samtools view -c -f 0x4 filtered.bam samtools view -c -f 0x100 filtered.bam samtools view -c -f 0x800 filtered.bam ``` Expected output: ```text 0 0 0 ``` If any count is nonzero, the corresponding bit was not in the exclusion mask. For example, a nonzero count for `0x100` means secondary alignments remain. Recheck the mask arithmetic and the command. Finally, compare flag statistics before and after filtering. `samtools flagstat` counts total records, primary records, secondary records, supplementary records, duplicates, mapped records, and paired records. The exact layout varies by version, but the categories are stable. ```bash samtools flagstat input.sorted.bam > input.flagstat.txt samtools flagstat filtered.bam > filtered.flagstat.txt diff -u input.flagstat.txt filtered.flagstat.txt | head -n 40 ``` Representative filtered output includes lines like these. ```text 742318 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 742318 + 0 mapped (100.00% : N/A) ``` The `+ 0` columns separate QC-passed and QC-failed records. If you excluded QC-failed reads, the second column should be zero for those categories. If you kept duplicates, the duplicate line will not be zero. ### Step 9: Inspect region and flag consistency Counts are necessary but not sufficient. Inspect a sample of the filtered alignments to confirm that the chromosomes and positions look like your target intervals. ```bash samtools view filtered.bam | awk 'BEGIN{OFS="\t"} {print $3,$4}' | sort -u | head ``` Representative output: ```text chr1 1001 chr2 50001 ``` The positions are 1-based leftmost alignment positions, not BED start coordinates. A read that starts at 1001 can overlap the first BED interval. A read that starts before 1001 can also overlap it if its alignment extends into the interval. That is why a simple position check is not a containment test. If you need reads fully contained inside intervals, `samtools view -L` is not enough. Use a region tool that can require full containment, or filter the output with a BED-aware tool and verify the overlap rule. For many analyses, overlap is the correct behavior because a read that partly covers the target still provides evidence at the target. ### Step 10: Move to downstream analysis A validated filtered BAM can go into variant calling, read counting, peak calling, visualization, or targeted extraction. Index the BAM before loading it into a genome browser. Keep the command line, the `samtools` version, the MAPQ cutoff, the flag mask, and the BED file with the output. Reproducibility depends on those details. For multithreaded SAM/BAM processing, elPrep is another tool designed for these files and for efficient memory use [1]. The same principles apply: know which records you are removing, validate the count, and keep the filter logic explicit. ## Practical Tips for Reliable BAM Filtering Write to a new file every time. Overwriting the input BAM makes it impossible to compare before and after counts or to recover records that you removed by mistake. Filter after duplicate marking if you plan to use the duplicate flag. The duplicate bit is not set by magic. A tool such as Picard MarkDuplicates or a equivalent marker must run first. Filtering for duplicates before marking them removes nothing. Keep a permissive BAM when the analysis may need it. A hard filter deletes records. If you might need multi-mapping reads, unmapped reads, or supplementary alignments later, keep an unfiltered or lightly filtered copy. Use `samtools view -c` for quick counts and `samtools flagstat` for category counts. The first counts alignment records. The second separates primary, secondary, supplementary, duplicate, mapped, and paired records. Use both when the distinction matters. Check chromosome naming before every region filter. The BAM header and the BED file must use the same names, including any `chr` prefix. A mismatch usually produces zero reads rather than a clear error. Record the exact flag mask in hexadecimal and in words. `-F 0xF04` is compact but opaque. A note such as "exclude unmapped, secondary, QC-fail, duplicate, supplementary" makes the command auditable months later. Do not assume that MAPQ 255 is a high-quality value. In many aligners, 255 means mapping quality is unavailable. A `-q 30` filter will keep MAPQ 255 records because 255 is greater than 30. If unavailable MAPQ values are unacceptable, inspect the distribution and handle those records explicitly. Watch memory use when you process many BAMs in parallel. Each `samtools view` process has its own buffers, and a high thread count multiplies memory use. Large SAM/BAM processing is memory-sensitive, which is why dedicated frameworks pay attention to memory management [1]. ## Common Mistakes and Limitations **Confusing MAPQ with base quality.** MAPQ is a placement confidence score. Base quality is a per-base error probability. Filtering a BAM with `-q 30` does not remove low-quality bases inside a read. It removes reads whose alignment placement is below the threshold. **Using one MAPQ cutoff across aligners.** BWA-MEM, Bowtie2, minimap2, and other aligners do not share a single calibrated MAPQ scale. A cutoff that works for one aligner can be too strict or too permissive for another. Inspect the distribution and read the aligner documentation. **Forgetting secondary and supplementary alignments.** A BAM can contain more than one alignment record for the same read. Secondary alignments have flag 0x100. Supplementary alignments have flag 0x800. If you exclude only unmapped reads, you can retain multiple records per read and inflate counts. Use `-F 0x900` to remove both classes when you need one primary alignment per read. **Forgetting duplicate and QC-fail flags.** Duplicates have flag 0x400. QC-failed reads have flag 0x200. These bits are set only if the upstream tool marked them. A duplicate filter before duplicate marking does nothing. A QC-fail filter before QC flagging does nothing. **Inverting `-f` and `-F`.** `-f 0x4` keeps unmapped reads. `-F 0x4` removes unmapped reads. The two commands differ by one letter and produce opposite results. Read the command aloud before you run it. **Counting alignments as read pairs.** `samtools view -c` counts alignment records. A paired-end fragment can yield two records. A BAM with 1,000,000 records may contain fewer than 1,000,000 fragments. Use `samtools flagstat` and, when needed, pair-aware counting tools to avoid overcounting. **Mixing BED and region-string coordinates.** BED is 0-based and half-open. A region string is 1-based and inclusive. A one-base shift is enough to change which reads overlap a boundary. Write down the convention before you build the BED file. **Assuming `-L` means full containment.** `samtools view -L` returns reads that overlap the intervals. A read can start outside the interval and still pass. If the analysis requires full containment, use a tool that enforces that rule or filter the output with explicit overlap logic. **Ignoring chromosome naming.** A BED file with `1` and a BAM header with `chr1` can produce zero reads without an error. Check the header with `samtools view -H` and confirm that every BED chromosome appears in the `@SQ` lines. **Using a stale or missing index.** Region queries on an indexed BAM need a current index. If you replace or rewrite the BAM, rebuild the index. A message such as `[E::idx_find_and_load] Could not retrieve index file for 'input.bam'` means the index is missing. Run `samtools index input.bam` after sorting. **Filtering before duplicate marking.** Duplicate marking compares read positions and orientations across the BAM. If you remove reads first, you change the duplicate-marking input and can change which reads are marked. Mark duplicates on the appropriate unfiltered or lightly filtered BAM, then filter on the duplicate flag. **Treating MAPQ 255 as a perfect score.** MAPQ 255 often means the aligner did not provide a mapping quality. A threshold of 30 keeps those records. If your analysis needs known MAPQ values, inspect and filter them explicitly. **Hard-filtering reads that downstream tools need.** Variant callers, [structural variant](/blog/guides/structural-variant) tools, and RNA-seq quantifiers may use multi-mapping reads, supplementary alignments, or unmapped reads in different ways. A hard filter removes those records permanently. Check the downstream tool's requirements before you set a strict cutoff. **Assuming a filtered BAM is analysis-ready.** A BAM that passes MAPQ, flag, and region filters can still have adapter contamination, overlapping mates, duplicate reads, or reference mismatches. Filtering is one step in a pipeline, not a substitute for library QC. **Memory and throughput limits.** Large BAMs can stress memory when many jobs run at once. Streaming commands help, but thread counts and parallel jobs still add up. SAM/BAM processing is memory-sensitive, and production pipelines need to budget for it [1]. **Biological bias from hard MAPQ filters.** Repeats, paralogs, and multi-copy genes often produce low MAPQ reads. Removing them can reduce ## Sources 1. [Comparing Ease of Programming in C++, Go, and Java for Implementing a Next-Generation Sequencing Tool.](https://pubmed.ncbi.nlm.nih.gov/31452597/)