🛠️ High-Throughput Alignment Workbench
SAMtools and BAMtools form the foundational standard for manipulating genomic, transcriptomic, and targeted sequence alignments in SAM, BAM, and CRAM formats. Use the interactive tools below to calculate bitwise flags, calculate sequencing coverage, and generate production-grade bash commands.
🚩 SAM Bitwise Flag Decoder & Command Generator
Interactive ToolEach mapped read in a SAM/BAM file has a 12-bit integer flag representing its paired-end state, orientation, and mapping validity. Check the properties below or type any flag integer (e.g. 99, 147, 163, 83, 1024) to decode it in real time.
samtools view -f 99 alignments.bam
📊 Coverage Depth & Lander-Waterman Estimator
Interactive CalculatorExperimental Parameters
Computed Metrics
samtools coverage -r chr1:1000000-2000000 input.bam
0. SAM Format & Header Specifications
SAM (Sequence Alignment/Map) files consist of optional header lines starting with @ followed by tab-delimited alignment records. Critical header tags include @HD (format version & sort order), @SQ (reference sequences & lengths), @RG (Read Groups linking reads to libraries/samples), and @PG (program provenance).
# Inspect SAM header without decompressing all records samtools view -H sample.bam # Extract and inspect specific Read Group tags (@RG) samtools view -H sample.bam | grep "^@RG"
1. SAM to BAM Compression & Multi-threaded Conversion
Converting human-readable SAM into binary compressed BGZF-blocked BAM reduces disk footprint by 75-80% and unlocks indexed random access.
# Convert SAM to BAM using 8 compression threads samtools view -@ 8 -bS raw_alignments.sam -o raw_alignments.bam # Validate BAM structural integrity and BGZF EOF block marker samtools quickcheck -v raw_alignments.bam
2. Coordinate vs Queryname Sorting
Downstream tools require specific sorting: Variant calling and genome browsers require coordinate sort (by chromosome and base position), while read pair fixers (fixmate) and count extractors (htseq-count) require queryname sort (by read identifier).
# Coordinate sort (default) allocating 2G RAM per thread samtools sort -@ 8 -m 2G -o sorted_coord.bam raw_alignments.bam # Queryname sort (grouped by read ID for pairing) samtools sort -n -@ 8 -m 2G -o sorted_name.bam raw_alignments.bam
3. BAM Indexing (Standard BAI vs CSI Large Chromosomes)
Indexing enables sub-second retrieval of reads spanning specific chromosomal coordinates without scanning gigabytes of data.
# Build standard BAM index (.bai for reference chromosomes < 512 Mbp) samtools index -@ 4 sorted_coord.bam # Build Coordinate-Sorted Index (.csi) for gigantic chromosomes (e.g. axolotl, wheat) samtools index -c -m 14 -@ 4 sorted_coord.bam
4. Alignment Metrics & Flag Statistics
Inspect overall alignment rate, mapped reads, properly paired proportions, and per-chromosome read counts.
# High-level flag metrics (mapped %, paired %, singletons) samtools flagstat -@ 4 sorted_coord.bam # Detailed per-chromosome read distribution (mapped vs unmapped counts) samtools idxstats sorted_coord.bam # Comprehensive alignment quality statistics (insert size distribution, error rate) samtools stats -@ 4 sorted_coord.bam > alignment_qc.stats
5. PCR & Optical Duplicate Marking (markdup)
PCR duplicates share identical 5' start coordinates and can introduce false allele frequency skews in variant calling. SAMtools provides a streamlined native markdup pipeline.
# Step 1: Queryname grouping with collate samtools collate -@ 8 -O -u sorted_coord.bam tmp-collate | \ # Step 2: Add ms and MC mate tags samtools fixmate -m -u - - | \ # Step 3: Coordinate sort samtools sort -@ 8 -m 2G -u - | \ # Step 4: Mark duplicates (add 0x400 flag) samtools markdup -@ 8 -s - marked_duplicates.bam
6. Targeted Regional Slicing & BED Filtering
Extract only alignments overlapping target genes, exons, or exome capture bait regions defined in a BED file.
# Slicing specific genomic interval samtools view -@ 4 -b sorted_coord.bam "chr7:140453136-140624564" -o braf_locus.bam # Slicing all regions listed in a BED capture panel samtools view -@ 4 -b -L exome_targets.bed sorted_coord.bam -o on_target.bam
7. CRAM Ultra-Compact Reference-Based Archiving
CRAM achieves 50-60% smaller file size than BAM by storing only differences relative to a reference FASTA file.
# Compress BAM to CRAM with reference genome samtools view -@ 8 -C -T reference.fasta -o archived.cram sorted_coord.bam # Index the CRAM archive samtools index archived.cram
8. Extracting Clean FASTQ Reads from BAM
Safely convert BAM back to paired-end FASTQ while properly routing singletons and discarding optical duplicates.
# Collate reads by queryname first, then extract cleanly samtools collate -u -O -@ 8 sorted_coord.bam | \ samtools fastq -@ 8 \ -1 read_1.fastq.gz \ -2 read_2.fastq.gz \ -0 /dev/null \ -s singleton.fastq.gz -
9. mpileup & BCFtools Variant Calling Pipeline
Generate genotype likelihoods across all reference positions and call SNVs and indels with the multiallelic caller.
# Stream mpileup into bcftools call for accurate SNV/INDEL calling bcftools mpileup -Ou -f reference.fasta -a FORMAT/AD,FORMAT/DP marked_duplicates.bam | \ bcftools call -mv -Oz -o raw_variants.vcf.gz # Index the resulting compressed VCF bcftools index -t raw_variants.vcf.gz
10. BAMtools Advanced Querying & Manipulation
Use BAMtools for JSON-based filtering, tag inspection, and split-by-chromosome workflows.
# Split BAM into individual chromosomes bamtools split -in sorted_coord.bam -reference # Count total alignments matching criteria bamtools count -in sorted_coord.bam -mapQuality ">=30"
🔀 Decision Architecture: Storage & Duplicate Strategy
Decision TreeCompress BAM to CRAM with
samtools view -C -T ref.fa. Saves ~55% storage. Always preserve exact reference FASTA checksums (MD5) to ensure lossless decoding.