WobblePlusAligner

Basic Usage

wobbleplusaligner.sh <query> <ref>
wobbleplusaligner.sh <query> <ref> <map>
wobbleplusaligner.sh <query> <ref> <map> <iterations> <simd>

WobblePlusAligner3 performs sequence alignment using a banded dynamic programming algorithm optimized for sequences with high identity and minimal indels. The tool uses only 2 arrays (prev[rLen+1] and curr[rLen+1]) avoiding full O(mn) matrix storage, implements traceback-free alignment calculation with 64-bit bit-packed encoding, and calculates bandwidth using Tools.mid(7, 1+max(qLen,rLen)/24, 24) with mismatch-based adjustment.

Parameters

Parameters are specified as positional arguments in the order shown below.

Examples

wobbleplusaligner.sh ATCGATCG ATCGATCG

wobbleplusaligner.sh query.fasta reference.fasta

wobbleplusaligner.sh query.fasta reference.fasta alignment_matrix.txt

wobbleplusaligner.sh query.fasta reference.fasta null 1000 simd

Algorithm Details

Core Algorithm

WobblePlusAligner3 implements a banded dynamic programming alignment algorithm using the decideBandwidth() method and RingBuffer for quality assessment. The algorithm encodes position and deletion information directly in 64-bit score values using bit fields (21-bit position, 21-bit deletion, 22-bit score), achieving O(n) space complexity through dual rolling arrays. Technical implementation:

Dynamic Bandwidth Calculation

The initial bandwidth is calculated using the formula: Tools.mid(7, 1+max(qLen,rLen)/24, 24), then refined by counting mismatches in sequence prefixes using decideBandwidth(). During alignment, bandwidth adapts using a RingBuffer(ringSize) where ringSize=(bandWidth0*5)/4 that tracks recent alignment quality. Poor local scores trigger bandwidth expansion up to ringSize*2 positions, while good scores maintain narrow bands to reduce computational work.

Band Center Tracking

The band center updates using: center = center + 1 + drift, where drift = Tools.mid(-1, maxPos-center, 2). This allows the band to follow the optimal alignment path while preventing drift beyond ±2 positions per row. The band boundaries are: bandStart = max(bandStart, center-bandWidth+quarterBand) and bandEnd = min(rLen, center+bandWidth+quarterBand).

Dual Array Matrix Processing

The algorithm allocates only two long arrays: prev[rLen+1] and curr[rLen+1], swapping roles each iteration. This eliminates the need for full O(n*m) matrix storage, reducing memory requirements to O(reference_length). Array initialization uses Arrays.fill(curr, BAD) where BAD = Long.MIN_VALUE/2 to prevent overflow during score calculations.

Long-Encoded Sequence Processing

Sequences are converted to long[] arrays using Factory.encodeLong() with N-base encoding: query N=15, reference N=31. The encoding formula ((q|r)>=15) detects N matches, which receive neutral scoring (N_SCORE=0L). This encoding enables vectorized operations and bitwise comparison during alignment matrix calculations using shared.SIMDAlign.alignBandVectorAndReturnMaxPos().

Scoring System

• Match: +1 point (1L << SCORE_SHIFT)
• Substitution: -1 point ((-1L) << SCORE_SHIFT)
• Insertion: -1 point ((-1L) << SCORE_SHIFT)
• Deletion: -1 point with position tracking
• N matches: 0 points (neutral scoring for ambiguous bases)

Bit-Field Score Encoding

Alignment information is encoded directly in 64-bit score values without requiring traceback matrices:

• POSITION_BITS (21): Reference start coordinate (bits 0-20, mask: 0x1FFFFF)
• DEL_BITS (21): Deletion count (bits 21-41, mask: 0x3FFFFFE00000)
• Score (22): Alignment score (bits 42-63, SCORE_SHIFT=42)

Score operations: MATCH = 1L << 42, SUB/INS/DEL = (-1L) << 42, DEL_INCREMENT = DEL + (1L << 21) for position tracking.

Vectorized Processing Mode

SIMD mode (enabled via command line flag) delegates band processing to shared.SIMDAlign.alignBandVectorAndReturnMaxPos(), which processes multiple alignment cells simultaneously using processor vector units. The vectorized code handles the main alignment loop (match/substitution/insertion calculations) while the Java code processes the deletion tail loop sequentially for correctness.

Quality-Based Bandwidth Expansion

A RingBuffer(ringSize) tracks recent maximum scores to detect local alignment degradation. When recentMissingScore = (oldScore + ringSize) - currentScore exceeds thresholds, scoreBonus = min(ringSize*2, recentMissingScore*2) expands bandwidth. Near sequence starts, additional expansion follows: bandWidth0 + max(10 + bandWidth0*8 - maxDrift*i, scoreBonus).

Performance Characteristics

• Time Complexity: O(n × w) where n is sequence length and w is average band width
• Space Complexity: O(r) where r is reference length
• Maximum Sequence Length: 2^21 positions (approximately 2 million bases)
• Optimal Use Case: High-identity sequences (>90%) with minimal indels

Identity Calculation

The algorithm extracts alignment statistics without traceback using bit field decoding:

• originPos: (maxScore & POSITION_MASK) extracts reference start
• endPos: maxPos provides reference end coordinate
• deletions: ((maxScore & DEL_MASK) >> POSITION_BITS) extracts deletion count
• rawScore: (maxScore >> SCORE_SHIFT) extracts alignment score

Identity calculation solves the system: (1) M+S+I=qLen, (2) M+S+D=refAlnLength, (3) Score=M-S-I-D, where:

• insertions: max(0, qLen + deletions - refAlnLength)
• matches: (rawScore + qLen + deletions) / 2.0
• substitutions: max(0, qLen - matches - insertions)
• identity: matches / (matches + substitutions + insertions + deletions)

Technical Specifications

Limitations

• Maximum reference sequence length: 2^21 - 1 positions (~2.1 million bases)
• Optimized for high-identity alignments (>80% identity)
• Performance may degrade for sequences with many large indels
• Does not support protein sequences (nucleotide sequences only)

Memory Usage

Memory requirements scale linearly with reference sequence length. Approximate memory usage:

• Base memory: ~16 bytes per reference position for score arrays
• Encoded sequences: ~8 bytes per position for long encoding
• Total: Approximately 24 bytes per reference position

Compatibility

• SIMD mode requires processor support for vectorized operations
• Visualization output compatible with visualizealignment.sh
• Supports both literal sequences and FASTA file input

Support

For questions and support:

• Email: bbushnell@lbl.gov
• Documentation: bbmap.org