2. Read Alignment & Mapping

BWA-MEM2 is the vector-accelerated (SSE, AVX2, AVX-512) implementation of the industry-benchmark BWA-MEM algorithm. It delivers 1.5–3x faster mapping of 70–250 bp paired-end reads to reference genomes while producing byte-identical alignments, serving as the gold standard for GATK germline and somatic workflows.

# 1. Index reference genome with BWA-MEM2
bwa-mem2 index -p GRCh38_bwa GRCh38.primary_assembly.genome.fa

# 2. Align paired-end FASTQ reads and stream directly to sorted BAM
REF="GRCh38_bwa"
RG="@RG\tID:FLOWCELL_L1\tSM:SAMPLE_01\tPL:ILLUMINA\tLB:LIB_DNA_1"

bwa-mem2 mem -t 16 -M -R "${RG}" "${REF}" \
  sample_clean_R1.fq.gz sample_clean_R2.fq.gz \
  | samtools sort -@ 4 -m 2G -o sample_aligned.sorted.bam -

# 3. Generate BAM index
samtools index -@ 4 sample_aligned.sorted.bam

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads of 50 up to several hundred base pairs to long reference genomes. It excels at local and end-to-end alignments, gapped alignments, and is the standard aligner in ENCODE ChIP-seq, Cut&Run, and ATAC-seq processing pipelines.

# 1. Build FM-index for reference genome
bowtie2-build --threads 8 GRCh38.fa hg38_bowtie2

# 2. Run highly sensitive end-to-end paired-end alignment
bowtie2 \
  -p 16 \
  --very-sensitive \
  -X 1000 \
  -x hg38_bowtie2 \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  2> sample_bowtie2_alignment.log \
  | samtools view -@ 4 -bS -F 4 -q 30 - \
  | samtools sort -@ 4 -o sample_bowtie2.sorted.bam -

# 3. Index BAM
samtools index sample_bowtie2.sorted.bam

Minimap2 is a versatile sequence alignment program that aligns DNA or mRNA sequences against a large reference database. It uses approximate minimizers for seeding followed by dynamic programming chaining, processing millions of noisy Oxford Nanopore or highly accurate PacBio HiFi reads at exceptional speed.

# 1. Oxford Nanopore Genomic DNA Alignment
minimap2 -ax map-ont -t 16 --MD -y \
  GRCh38.fa ont_reads.fastq.gz \
  | samtools sort -@ 4 -m 2G -o ont_aligned.sorted.bam -

# 2. PacBio HiFi Genomic DNA Alignment
minimap2 -ax map-hifi -t 16 --MD \
  GRCh38.fa pacbio_hifi.fastq.gz \
  | samtools sort -@ 4 -m 2G -o pacbio_aligned.sorted.bam -

# 3. Direct RNA / full-length cDNA spliced alignment
minimap2 -ax splice -uf -k14 -t 16 \
  GRCh38.fa ont_cdna.fastq.gz \
  | samtools sort -@ 4 -o ont_spliced_rna.sorted.bam -

samtools index ont_aligned.sorted.bam

Spliced Transcripts Alignment to a Reference (STAR) aligns non-contiguous RNA-seq reads across canonical and non-canonical splice junctions using sequential maximum mappable seed searches. It is the core spliced aligner used by the ENCODE Consortium, GTEx, and the nf-core/rnaseq standard pipeline.

# 1. Generate STAR Genome Index with GTF Annotation
STAR --runMode genomeGenerate --runThreadN 16 \
  --genomeDir ./star_index/ \
  --genomeFastaFiles GRCh38.primary_assembly.fa \
  --sjdbGTFfile gencode.v46.annotation.gtf \
  --sjdbOverhang 149

# 2. Run 2-Pass Spliced Alignment with Direct Sorted BAM Output
STAR --runThreadN 16 \
  --genomeDir ./star_index/ \
  --readFilesIn sample_R1.fq.gz sample_R2.fq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --twopassMode Basic \
  --outFileNamePrefix ./star_out/sample_ \
  --quantMode GeneCounts

# Index the sorted alignment
samtools index star_out/sample_Aligned.sortedByCoord.out.bam

Hierarchical Indexing for Spliced Alignment of Transcripts (HISAT2) uses a Hierarchical Graph FM index that represents both the reference genome and millions of known genomic variants and splice junctions. It operates with a tiny RAM footprint (<8 GB), making it suitable for memory-constrained workstations.

# 1. Extract splice sites and exons from GTF
hisat2_extract_splice_sites.py gencode.v46.gtf > splicesites.tsv
hisat2_extract_exons.py gencode.v46.gtf > exons.tsv

# 2. Build graph-based reference index
hisat2-build -p 16 --ss splicesites.tsv --exon exons.tsv GRCh38.fa hisat2_hg38

# 3. Execute paired-end spliced alignment with strand orientation
hisat2 -p 16 \
  --rna-strandness RF \
  --known-splicesite-infile splicesites.tsv \
  -x hisat2_hg38 \
  -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
  --summary-file hisat2_summary.txt \
  | samtools sort -@ 4 -o hisat2_aligned.sorted.bam -

samtools index hisat2_aligned.sorted.bam

pbmm2 is Pacific Biosciences’ official Minimap2 wrapper tailored for PacBio long-read data. It provides verified preset parameter weighting for HiFi, CCS, and Kinnex transcript reads while preserving internal BAM read groups and kinematic kinetic tags required for downstream methylation calling.

# 1. Align unaligned PacBio BAM directly to reference FASTA
pbmm2 align \
  GRCh38.fa \
  m84011_230905_sample.hifi_reads.bam \
  sample_pbmm2_aligned.sorted.bam \
  --preset HIFI \
  --sort \
  -j 16 \
  --log-level INFO

# 2. View alignment summary statistics
samtools flagstat sample_pbmm2_aligned.sorted.bam

Winnowmap2 addresses the fundamental challenge of mapping long reads to repetitive genomic regions, centromeric arrays, and segmental duplications. By calculating k-mer frequencies and downweighting repetitive minimizers, it prevents ambiguous anchors and drastically improves structural variant and variant calling in complex loci.

# 1. Pre-calculate repetitive 15-mers using meryl
meryl count k=15 output merylDB GRCh38.fa
meryl print greater-than distinct=0.9998 merylDB > repetitive_k15.txt

# 2. Run Winnowmap2 with repeat-downweighting
winnowmap -W repetitive_k15.txt \
  -ax map-ont -t 16 --MD \
  GRCh38.fa \
  ont_ultralong_reads.fastq.gz \
  | samtools sort -@ 4 -m 2G -o winnowmap_aligned.sorted.bam -

samtools index winnowmap_aligned.sorted.bam