3. BAM/CRAM Processing & General Utilities

SAMtools provides the essential command-line infrastructure for reading, writing, sorting, indexing, and calculating coverage statistics on high-throughput sequencing alignments in SAM, BAM, and CRAM formats.

# 1. Multi-threaded coordinate sort and indexing
samtools sort -@ 8 -m 2G -o sample.sorted.bam sample.raw.bam
samtools index -@ 4 sample.sorted.bam

# 2. Extract alignment flag statistics
samtools flagstat -@ 4 sample.sorted.bam > sample.flagstat.txt
samtools idxstats sample.sorted.bam > sample.idxstats.txt

# 3. Filter for high-quality uniquely mapped paired reads (MAPQ >= 30)
samtools view -@ 8 -b -f 2 -F 3844 -q 30 \
  sample.sorted.bam chr1 chr2 chr3 \
  -o sample.filtered.bam

# 4. Convert BAM to space-saving CRAM format with reference
samtools view -@ 8 -C -T GRCh38.fa sample.filtered.bam -o sample.cram
samtools index sample.cram

Picard is a set of command-line tools for manipulating high-throughput sequencing data in BAM/CRAM formats. Its MarkDuplicates module is the standard GATK best practice for flagging PCR and optical duplicate molecules.

# 1. Mark optical & PCR duplicate reads with Picard
picard MarkDuplicates \
  I=sample.sorted.bam \
  O=sample.dedup.bam \
  M=sample.duplicate_metrics.txt \
  OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 \
  CREATE_INDEX=true \
  VALIDATION_STRINGENCY=SILENT

# 2. Collect Insert Size Metrics for Paired-End Library QC
picard CollectInsertSizeMetrics \
  I=sample.dedup.bam \
  O=sample.insert_size_metrics.txt \
  H=sample.insert_size_histogram.pdf \
  M=0.5

BEDtools is an extensive suite of utilities for genomic interval manipulation. It enables comparing, intersecting, merging, and extracting sequence features across BED, GFF/GTF, VCF, and BAM files with single-base exactness.

# 1. Intersect BAM alignments with whole-exome capture targets
bedtools intersect -abam sample.dedup.bam -b exome_targets.bed > sample.on_target.bam

# 2. Merge overlapping peak intervals in ChIP/ATAC datasets
bedtools merge -i sample_peaks.bed -c 4,5 -o count,mean > merged_peaks.bed

# 3. Compute per-base physical sequencing coverage across custom intervals
bedtools coverage -a candidate_genes.bed -b sample.dedup.bam > gene_coverage.txt

# 4. Extract flanking 500bp promoter regions upstream of TSS
bedtools flank -i tss.bed -g hg38.chrom.sizes -l 500 -r 0 -s > promoters_500bp.bed

Sambamba is a high-performance alternative to SAMtools and Picard written in D. It parallelizes computationally intensive operations like BAM sorting, indexing, and duplicate marking, reducing runtime by up to 5x on multicore HPC systems.

# 1. High-speed multi-threaded coordinate sort
sambamba sort -t 16 -m 8GB -o sample.sambamba_sorted.bam sample.raw.bam

# 2. Fast parallel duplicate marking
sambamba markdup -t 16 --overflow-list-size 600000 \
  sample.sambamba_sorted.bam sample.sambamba_dedup.bam

# 3. Rapid slice querying and filtering
sambamba view -t 8 -f bam -F "mapping_quality >= 30 and not (unmapped or duplicate)" \
  sample.sambamba_dedup.bam -o sample.clean.bam

deepTools is a user-friendly suite of Python tools developed for analyzing and visualizing deep sequencing data (ChIP-seq, ATAC-seq, RNA-seq). It creates genome-wide 1x-normalized bigWig files, correlation heatmaps, and publication-quality TSS meta-gene profiles.

# 1. Generate RPGC (1x coverage) normalized bigWig file
bamCoverage -b sample.dedup.bam -o sample.rpgc.bw \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --binSize 10 \
  --smoothLength 30 \
  --numberOfProcessors 12

# 2. Compute signal matrix around Transcription Start Sites (TSS)
computeMatrix reference-point \
  --referencePoint TSS \
  -b 2000 -a 2000 \
  -R gencode.v46.protein_coding_genes.bed \
  -S sample.rpgc.bw \
  -o matrix_tss.gz \
  --numberOfProcessors 12

# 3. Render publication-ready profile and heatmap
plotHeatmap -m matrix_tss.gz \
  -out tss_enrichment_heatmap.pdf \
  --colorMap viridis \
  --whatToShow 'plot, heatmap and colorbar'

mosdepth is a blazing fast BAM/CRAM depth calculation tool written in Nim. It computes per-base, windowed, and region-level mean depth across entire mammalian genomes in minutes, outputting threshold distributions (e.g. % bases >= 20x, 30x) required for clinical genomics validation.

# 1. Calculate target coverage and clinical thresholds on Exome BAM
mosdepth --threads 8 \
  --by exome_targets.bed \
  --thresholds 10,20,30,50,100 \
  --fast-mode \
  sample_exome_depth \
  sample.dedup.bam

# 2. View summary table directly
cat sample_exome_depth.mosdepth.summary.txt

# 3. Whole Genome Sequencing 500bp windowed coverage
mosdepth --threads 8 --by 500 sample_wgs_500bp sample_wgs.cram