WFA-JS/pkg/wfa.go

package wfa

func WFAlign(s1 string, s2 string, penalties Penalty, doCIGAR bool) Result {
	n := len(s1)
	m := len(s2)
	A_k := m - n
	A_offset := m
	score := 0
	estimatedScore := (max(n, m) * max(penalties.M, penalties.X, penalties.O, penalties.E)) / 4
	M := NewWavefrontComponent(estimatedScore)
	M.SetLoHi(0, 0, 0)
	M.SetVal(0, 0, 0)
	M.SetTraceback(0, 0, End)
	I := NewWavefrontComponent(estimatedScore)
	D := NewWavefrontComponent(estimatedScore)

	for {
		WFExtend(M, s1, n, s2, m, score)
		ok, val := M.GetVal(score, A_k)
		if ok && val >= A_offset {
			break
		}
		score = score + 1
		WFNext(M, I, D, score, penalties)
	}

	CIGAR := ""
	if doCIGAR {
		CIGAR = WFBacktrace(M, I, D, score, penalties, A_k, s1, s2)
	}

	return Result{
		Score: score,
		CIGAR: CIGAR,
	}
}

func WFExtend(M WavefrontComponent, s1 string, n int, s2 string, m int, score int) {
	_, lo, hi := M.GetLoHi(score)
	for k := lo; k <= hi; k++ {
		// v = M[score][k] - k
		// h = M[score][k]
		ok, h := M.GetVal(score, k)
		v := h - k

		// exit early if v or h are invalid
		if !ok {
			continue
		}
		for v < n && h < m && s1[v] == s2[h] {
			_, val := M.GetVal(score, k)
			M.SetVal(score, k, val+1)
			v++
			h++
		}
	}
}

func WFNext(M WavefrontComponent, I WavefrontComponent, D WavefrontComponent, score int, penalties Penalty) {
	// get this score's lo, hi
	lo, hi := NextLoHi(M, I, D, score, penalties)

	for k := lo; k <= hi; k++ {
		NextI(M, I, score, k, penalties)
		NextD(M, D, score, k, penalties)
		NextM(M, I, D, score, k, penalties)
	}
}

func WFBacktrace(M WavefrontComponent, I WavefrontComponent, D WavefrontComponent, score int, penalties Penalty, A_k int, s1 string, s2 string) string {
	traceback_CIGAR := []string{"I", "I", "D", "D", "X", "", "", ""}
	x := penalties.X
	o := penalties.O
	e := penalties.E
	CIGAR_rev := ""
	tb_s := score
	tb_k := A_k
	_, current_traceback := M.GetTraceback(tb_s, tb_k)
	done := false

	for !done {
		CIGAR_rev = CIGAR_rev + traceback_CIGAR[current_traceback]
		switch current_traceback {
		case OpenIns:
			tb_s = tb_s - o - e
			tb_k = tb_k - 1
			_, current_traceback = M.GetTraceback(tb_s, tb_k)
		case ExtdIns:
			tb_s = tb_s - e
			tb_k = tb_k - 1
			_, current_traceback = I.GetTraceback(tb_s, tb_k)
		case OpenDel:
			tb_s = tb_s - o - e
			tb_k = tb_k + 1
			_, current_traceback = M.GetTraceback(tb_s, tb_k)
		case ExtdDel:
			tb_s = tb_s - e
			tb_k = tb_k + 1
			_, current_traceback = D.GetTraceback(tb_s, tb_k)
		case Sub:
			tb_s = tb_s - x
			// tb_k = tb_k;
			_, current_traceback = M.GetTraceback(tb_s, tb_k)
		case Ins:
			// tb_s = tb_s;
			// tb_k = tb_k;
			_, current_traceback = I.GetTraceback(tb_s, tb_k)
		case Del:
			// tb_s = tb_s;
			// tb_k = tb_k;
			_, current_traceback = D.GetTraceback(tb_s, tb_k)
		case End:
			done = true
		}
	}

	CIGAR_part := Reverse(CIGAR_rev)
	c := 0
	i := 0
	j := 0
	for i < len(s1) && j < len(s2) {
		if s1[i] == s2[j] {
			//CIGAR_part.splice(c, 0, "M")
			CIGAR_part = Splice(CIGAR_part, 'M', c)
			c++
			i++
			j++
		} else if CIGAR_part[c] == 'X' {
			c++
			i++
			j++
		} else if CIGAR_part[c] == 'I' {
			c++
			j++
		} else if CIGAR_part[c] == 'D' {
			c++
			i++
		}
	}

	return CIGAR_part
}