sequence-algorithms/figures/part2/lcsq.lua

function lcsq_matrix(seq1, seq2)
    local gap_penalty = 0
    local match_score = 1
    local n1 = string.len(seq1)
    local n2 = string.len(seq2)
    -- Create a n1 x n2 matrix
    local matrix = {}
    for i=0,n1 do
        matrix[i] = {}
        for j=0,n2 do
            matrix[i][j] = 0
        end
    end
    -- Fill the rest of the matrix
    local match, delete, insert
    for i=1,n1 do
        for j=1,n2 do
            if string.sub(seq1, i, i) == string.sub(seq2, j, j) then
                match = matrix[i-1][j-1] + match_score
            else
                match = matrix[i-1][j-1]
            end
            gap1 = matrix[i-1][j] + gap_penalty
            gap2 = matrix[i][j-1] + gap_penalty
            matrix[i][j] = math.max(match, gap1, gap2)
        end
    end
    return matrix
end

local function has_value (tab, val)
    for index, value in ipairs(tab) do
        if value == val then
            return true
        end
    end

    return false
end

function repr_matrix(matrix)
  repr = ""
  for i=0,#matrix do
    for j=0,#matrix[i] do
      repr = repr .. matrix[i][j] .. " "
    end
    repr = repr .. "\n"
  end
  return repr
end


function draw_lcsq_matrix_graph(seq1, seq2, matrix)
    local tikz_code = ""
    function coordinate(i, j)
        return i .. "_" .. j
    end
    local steps = {
        {-1, -1},
        {0, -1},
        {-1, 0},
    }

    local n1 = string.len(seq1)
    local n2 = string.len(seq2)
    local path = {}
    local i = n1
    local j = n2
    while i >= 0 and j >= 0 do
        path[#path+1] = coordinate(i, j)
        local max = matrix[i][j]
        local max_step = steps[1]
        for index, step in ipairs(steps) do
            local k = i + step[1]
            local l = j + step[2]
            if k >= 0 and l >= 0 and matrix[k][l] > max then
                max_step = step
                max = matrix[k][l]
            end
        end
        i = i + max_step[1]
        j = j + max_step[2]
    end
    -- Draw the matrix as tikz node with matrix value
    for i=0,n1 do
        for j=0,n2 do
            local options = ""
            if has_value(path, coordinate(i, j)) then

                options = "[fill=gray, draw, minimum size=1]"
            end
            tikz_code = tikz_code .. "\\node" .. options .. " (" .. coordinate(i, j) .. ") at (" .. i .. ", " .. -j .. ")" .. " {" .. matrix[i][j] .. "};"
        end
    end
    -- Add nucleotide labels
    for i=1,n1 do
        local nt = string.sub(seq1, i, i)
        tikz_code = tikz_code .. "\\node at (".. i .. "," .. 1 .. ")" .. "{$" .. nt .."$};"
    end
    for i=1,n2 do
        local nt = string.sub(seq2, i, i)
        tikz_code = tikz_code .. "\\node at (" .. -1 .. ", " .. -i .. ")" .. "{$ ".. nt .."$};"
    end
    -- For seq2
       for i=0,n1 do
        for j=0,n2 do
            local max = 0
            for index, step in ipairs(steps) do
                local k = i + step[1]
                local l = j + step[2]
                if k >= 0 and l >= 0 and matrix[k][l] > max then
                    max = matrix[k][l]
                end
            end
            local highlighted = false
            for index, step in ipairs(steps) do
                local k = i + step[1]
                local l = j + step[2]
                if k >= 0 and l >= 0 and matrix[k][l] == max then
                    tikz_code = tikz_code .. "\\draw[->] (" .. coordinate(i, j) .. ")" .. " -- " .. "(" .. coordinate (k, l) .. ");"
                end
            end
        end
    end
    return tikz_code
end

function draw_lcsq_matrix(seq1, seq2)
    -- print(string.format(" Path: %s -> %s", seq1, seq2))
    local matrix = lcsq_matrix(seq1, seq2)
    local n1 = string.len(seq1)
    local n2 = string.len(seq2)
    local repr = ""
    -- Draw the matrix as tikz nodes
    for i=0,n1-1 do
        for j=0,n2-1 do
            repr = repr .. " " .. string.format("\\node[draw, minimum width=1cm, minimum height=1cm] at (%d, -%d) {};", i, j, matrix[i][j])
        end
    end
    -- Draw the sequence labels
    for i=1,n1 do
        repr = repr .. " " .. string.format("\\node at (%d, -%d) {%s};", i-1, -1, string.sub(seq1, i, i))
    end
    for i=1,n2 do
      repr = repr .. " " .. string.format("\\node at (%d, -%d) {%s};", -1, i-1, string.sub(seq2, i, i))
    end
    -- Add a path from the bottom right corner to the top left corner, following the minimum of the three possible moves at each step
    local i, j, value, previous_value
    i = n1-1
    j = n2-1
    repr = repr .. " " string.format("\\draw[-,line width=2, gray] (%d, -%d) --", i, j)
    while i > 0 and j > 0 do
        value = math.min(matrix[i-1][j-1], matrix[i-1][j], matrix[i][j-1])
        if value == matrix[i-1][j-1] then
            i = i - 1
            j = j - 1
        elseif value == matrix[i-1][j] then
            i = i - 1
        else
            j = j - 1
        end
        repr = repr .. " " .. string.format(" (%d, -%d) -- ", i, j)
    end
    repr = repr .. " " .. string.format("(0, 0) -- (-1, 1);", i, j)
    return repr
end

function main()
  local seq1 = "ATCTGAT"
  local seq2 = "TGCATA"
  local matrix = lcsq_matrix(seq1, seq2)
  print(draw_lcsq_matrix_graph(seq1, seq2, matrix))
end

-- main()

return {
  lcsq_matrix=lcsq_matrix,
  draw_lcsq_matrix_graph=draw_lcsq_matrix_graph,
  draw_lcsq_matrix=draw_lcsq_matrix
}