Wavefront alignment algorithm (WFA) in Golang

This golang package implements Wavefront alignment algorithm (WFA), not BiWFA (maybe in the future).

Implemented features:

• Distance metrics: gap-affine.
• Alignment types: global, semi-global.
• Heuristics: wf-adaptive.

Executable binaries for common operating systems are available to download, for benchmarking with lots of sequences, or for fast aligning two sequences in the command line.

$ wfa-go -g "Bioinformatics helps Biology" "We learn bioinformatics to help biologists"
query   ---------Bioinformatics ---helps Biology---
                  ||||||||||||||   |||| | |||||
target  We learn bioinformatics to help- biologists
cigar   9I1X14M3I4M1D1M1X5M1X3I

align-score : 32
match-region: q[2, 27]/28 vs t[11, 38]/42
align-length: 29, matches: 24 (82.76%), gaps: 4, gap regions: 2

Table of Contents

• History
• Details
• Examples
• Usages
• CLI
• Benchmark
• Reference

History

1. I need a fast DNA alignment package in Golang for my project.
2. WFA "looks easy" to implement as it does not heavily reply on SIMD intrinsics, though there are some other algorithms performing well in a benchmark.
3. I'm not familiar with C++, and I found it difficult to understand the official code. I only found one (in Rust) 3rd party implementation.
4. After reading the WFA paper, I thought the algorithm is easy, so I implemented WFA from the scratch.
5. Later I found it not easy, there were so many details. After reading the next step in the rust version, I realised it's because I don't know gap-affine penalties.
6. The backtrace step is the most difficult part. In v0.1.0, I used 3 extra components to store the source offsets for each cell in I, D, M components. But it needs more memory. And the speed is also not ideal, 1/20 of official version.
   • Besides, I checked the bases again when tracing back matches. WFA did this too, but WFA2 did not.
7. Next, aftering reading thofficial implementation, I rewrote the whole project, using similar backtrace workfow with WFA2. The speed increased to 1/10 of the official version.
8. C++ is wild, it even support accessing list/array elements with negative indexes (the diagonal k).

Details

• A WFA component is saved as a list of WaveFront []*WaveFront.

  • Score s is the index.
  • If the value is nil, it means the score does not exsit.

• A WaveFront is saved as a list of offsets []uint32. We preset the list with a big length (2048) to avoid frequent append operations.

  • Diagonal k is the index. To support negative k values, we use ZigZag encoding:

    index: 0,  1,  2,  3,  4,  5,  6
    k:     0, -1,  1, -2,  2, -3,  3
    

    Value 0 means the k does not exist.

  • Offsets are saved with uint32 integers, with the lower 3 bits for saving 6 possible paths which are used for backtrace.

    wfaInsertOpen uint32 = iota + 1
    wfaInsertExt
    wfaDeleteOpen
    wfaDeleteExt
    wfaMismatch
    wfaMatch
    

• Maximum sequence length: 512 Mb, 536,870,911 (1<<(32-3) - 1).

• All objects are saved in object pools for computation efficiency. Just don't forget to recycle them.

Examples

Each WFA component can be visualized as a table. A table cell contains the alignment type symbol and the score.

⊕    Unknown,
⟼    Gap open (Insertion)
🠦    Gap extension (Insertion)
↧    Gap open (Deletion)
🠧    Gap extension (Deletion)
⬂    Mismatch
⬊    Match

Global alignment

CIGAR:  1M2X2M1X4M
query   ACCATACTCG
        |  || ||||
target  AGGATGCTCG

align-score : 12
align-region: q[1, 10] vs t[1, 10]
align-length: 10, matches: 7 (70.00%), gaps: 0, gapRegions: 0

Semi-global alignment

CIGAR:  1I1M1X1M1X1M1I4M1I
query   -ACGAT-CTCG-
         | | | ||||
target  CAGGCTCCTCGG

align-score : 16
align-region: q[1, 9] vs t[2, 11]
align-length: 10, matches: 7 (70.00%), gaps: 1, gapRegions: 1

Usages

import "github.com/shenwei356/wfa"

// ------------------[ initialization ]------------------

// aligner
algn := wfa.New(
    &wfa.Penalties{
        Mismatch: 4,
        GapOpen:  6,
        GapExt:   2,
    },
    &wfa.Options{
        GlobalAlignment: false,
    })

// set adaptive reduction parameters
algn.AdaptiveReduction(&wfa.AdaptiveReductionOption{
    MinWFLen:    10,
    MaxDistDiff: 50,
    CutoffStep:  1,
})

// ------------------[ align one pair of seqs ]------------------

q := []byte("ACCATACTCG")
t := []byte("AGGATGCTCG")

// align
result, err := algn.Align(q, t)
checkErr(err)

// score table of M
algn.Plot(&q, &t, os.Stdout, algn.M, true)

if outputAlignment {
    fmt.Println()
    fmt.Printf("CIGAR:  %s\n", result.CIGAR())

    Q, A, T := result.AlignmentText(&q, &t)
    fmt.Printf("query   %s\n", *Q)
    fmt.Printf("        %s\n", *A)
    fmt.Printf("target  %s\n", *T)

    fmt.Println()
    fmt.Printf("align-score : %d\n", result.Score)
    fmt.Printf("align-region: q[%d, %d] vs t[%d, %d]\n",
        result.QBegin, result.QEnd, result.TBegin, result.TEnd)
    fmt.Printf("align-length: %d, matches: %d (%.2f%%), gaps: %d, gapRegions: %d\n",
        result.AlignLen, result.Matches, float64(result.Matches)/float64(result.AlignLen)*100,
        result.Gaps, result.GapRegions)
    fmt.Println()

    
    wfa.RecycleAlignmentText(Q, A, T) // !! important, recycle objects
}

wfa.RecycleAlignmentResult(result) // !! important, recycle objects

// ------------------[ clean ]------------------

wfa.RecycleAligner(algn) // !! important, recycle objects

CLI

A CLI (download binaries) is available to align two sequences from either positional arguments or an input file (format).

$ wfa-go AGCTAGTGTCAATGGCTACTTTTCAGGTCCT AACTAAGTGTCGGTGGCTACTATATATCAGGTCCT
query   AGCTA-GTGTCAATGGCTACT---TTTCAGGTCCT
        | ||| |||||  ||||||||   | |||||||||
target  AACTAAGTGTCGGTGGCTACTATATATCAGGTCCT
cigar   1M1X3M1I5M2X8M3I1M1X9M

align-score : 36
match-region: q[1, 31]/31 vs t[1, 35]/35
align-length: 35, matches: 27 (77.14%), gaps: 4, gap regions: 2

$ wfa-go -i seqs.txt
query   A-TTGGAAAATAGGATTGGGGTTTGTTTATATTTGGGTTGAGGGATGTCCCACCTTCGTCGTCCTTACGTTTCCGGAAGGGAGTGGTTAGCTCGAAGCCCA
          |||||||||||||| ||||||||||||||||||||||||||||||||||||||| ||||||||||||||||||||||||||||||| ||||||||||||
target  GATTGGAAAATAGGAT-GGGGTTTGTTTATATTTGGGTTGAGGGATGTCCCACCTT-GTCGTCCTTACGTTTCCGGAAGGGAGTGGTT-GCTCGAAGCCCA
cigar   1X1I14M1D39M1D31M1D12M

align-score : 36
match-region: q[2, 100]/100 vs t[3, 98]/98
align-length: 99, matches: 96 (96.97%), gaps: 3, gap regions: 3

WFA alignment in Golang

 Author: Wei Shen <[email protected]>
   Code: https://github.com/shenwei356/wfa
Version: v0.2.0

Input file format:
  see https://github.com/smarco/WFA-paper?tab=readme-ov-file#41-introduction-to-benchmarking-wfa-simple-tests
  Example:
  >ATTGGAAAATAGGATTGGGGTTTGTTTATATTTGGGTTGAGGGATGTCCCACCTTCGTCGTCCTTACGTTTCCGGAAGGGAGTGGTTAGCTCGAAGCCCA
  <GATTGGAAAATAGGATGGGGTTTGTTTATATTTGGGTTGAGGGATGTCCCACCTTGTCGTCCTTACGTTTCCGGAAGGGAGTGGTTGCTCGAAGCCCA
  >CCGTAGAGTTAGACACTCGACCGTGGTGAATCCGCGACCACCGCTTTGACGGGCGCTCTACGGTATCCCGCGATTTGTGTACGTGAAGCAGTGATTAAAC
  <CCTAGAGTTAGACACTCGACCGTGGTGAATCCGCGATCTACCGCTTTGACGGGCGCTCTACGGTATCCCGCGATTTGTGTACGTGAAGCGAGTGATTAAAC

Usage:
  1. Align two sequences from the positional arguments.

        wfa-go [options] <query seq> <target seq>

  2. Align sequence pairs from the input file (described above).

        wfa-go [options] -i input.txt

Options/Flags:
  -N    do not output alignment (for benchmark)
  -a    do not use adaptive reduction
  -g    do not use global alignment
  -h    print help message
  -i string
        input file.
  -m    mem pprof. go tool pprof -http=:8080 mem.pprof
  -p    cpu pprof. go tool pprof -http=:8080 cpu.pprof

Benchmark

Generate datasets with WFA (e634175) or WFA2-lib (v2.3.5):

./bin/generate_dataset -n 100000 -l 1000 -e 0.05 -o l1000-e0.05.seq
./bin/generate_dataset -n 100000 -l 1000 -e 0.10 -o l1000-e0.10.seq
./bin/generate_dataset -n 100000 -l 1000 -e 0.20 -o l1000-e0.20.seq

./bin/generate_dataset -n 500 -l 50000 -e 0.05 -o l50000-e0.05.seq
./bin/generate_dataset -n 500 -l 50000 -e 0.10 -o l50000-e0.10.seq
./bin/generate_dataset -n 500 -l 50000 -e 0.20 -o l50000-e0.20.seq

Commands (all tools use wfa-adaptive heuristic 10,50,1):

# memusg: https://github.com/shenwei356/memusg

# WFA
memusg -t -s "./bin/align_benchmark -i l1000-e0.05.seq -a gap-affine-wfa-adaptive"

# WFA2-lib
memusg -t -s "./bin/align_benchmark -i l1000-e0.05.seq -a gap-affine-wfa --wfa-heuristic wfa-adaptive --wfa-heuristic-parameters 10,50,1"

# WFA-go (this package, binary files are availabe in the release page)
# global alignment && do not output results
memusg -t -s "wfa-go -N -i l1000-e0.05.seq"

csvtk csv2md -t benchmark.tsv -a c,c,c,l,r,r

Results:

Run in a laptop PC, with single-thread.

Reference

• Santiago Marco-Sola, Juan Carlos Moure, Miquel Moreto, Antonio Espinosa. "Fast gap-affine pairwise alignment using the wavefront algorithm." Bioinformatics, 2020.
• Santiago Marco-Sola, Jordan M Eizenga, Andrea Guarracino, Benedict Paten, Erik Garrison, Miquel Moreto. "Optimal gap-affine alignment in O(s) space". Bioinformatics, 2023.
• https://github.com/smarco/WFA-paper/
• https://github.com/smarco/WFA2-lib
• https://github.com/rchikhi/rust-alignbench

Support

Please open an issue to report bugs, propose new functions or ask for help.

License

MIT License