(The entire project lives here.)

Problem definition

Suppose we are given three genomic strings drawn from the alphabet of 20 amino acids:

ACGH

CFG

EAC

In multiple sequence alignment problem we want to align them optimally in order to access their evolutionary similarity. Such an alignment may look like the following:

-AC-GH

--CFG-

EAC---

Solution

Often, the problem is cast as a pathfinding problem in $k$-dimensional lattice, where $k$ is the number of strings. Here is my program, it compares A* with Dijkstra's algorithms:

MultipleSequenceAlignmentInstance.java

package net.coderodde.bio.msa;



import java.util.ArrayList;

import java.util.Collections;

import java.util.HashMap;

import java.util.HashSet;

import java.util.List;

import java.util.Map;

import java.util.Objects;

import java.util.PriorityQueue;

import java.util.Queue;

import java.util.Set;



public final class MultipleSequenceAlignmentInstance {



    /**

     * The penalty for pairs in which there is one valid character and one gap.

     */

    private final int gapPenalty;



    /**

     * The character pairs cost matrix.

     */

    private final CostMatrix<Integer> costMatrix;



    /**

     * The array of sequences to be aligned.

     */

    private final String sequenceArray;



    // A small speed optimization:

    private final char column;



    public MultipleSequenceAlignmentInstance(CostMatrix<Integer> costMatrix,

                                             int gapPenalty,

                                             String... sequenceArray) {

        this.costMatrix = Objects.requireNonNull(costMatrix,

                                                 "Cost matrix is null");

        this.gapPenalty = gapPenalty;

        this.sequenceArray = sequenceArray.clone();

        this.column = new char[sequenceArray.length];



        for (int i = 0; i != sequenceArray.length; ++i) {

            this.sequenceArray[i] = 

                    checkIsValidGenomicSequence(sequenceArray[i]);

        }

    }



    public Alignment align() {

        long start = System.currentTimeMillis();

        HeuristicFunction hf = new HeuristicFunctionComputer()

                .computeHeuristicFunction(this);

        long end = System.currentTimeMillis();



        System.out.println("Computed heuristic function in " + (end - start) +

                           " milliseconds.");



        Queue<LatticeNodeHolder> open = new PriorityQueue<>();

        Map<LatticeNode, Integer> distance = new HashMap<>();

        Map<LatticeNode, LatticeNode> parents = new HashMap<>();



        LatticeNode sourceNode = getSourceNode();

        LatticeNode targetNode = getTargetNode();



        open.add(new LatticeNodeHolder(sourceNode, 0));

        distance.put(sourceNode, 0);

        parents.put(sourceNode, null);



        Set<LatticeNode> closed = new HashSet<>();



        while (true) {

            LatticeNode currentNode = open.remove().getNode();



            if (currentNode.equals(targetNode)) {

                return tracebackPath(parents, distance.get(currentNode));

            }



            if (closed.contains(currentNode)) {

                continue;

            }



            closed.add(currentNode);



            for (LatticeNode childNode : currentNode.getChildren()) {

                if (closed.contains(childNode)) {

                    continue;

                }



                int tentativeCost = distance.get(currentNode) +

                                    getWeight(currentNode, childNode);

                Integer currentCost = distance.get(childNode);



                if (currentCost == null || currentCost > tentativeCost) {

                    open.add(new LatticeNodeHolder(childNode, 

                                                   tentativeCost + 

                                                           hf.get(childNode)));

                    distance.put(childNode, tentativeCost);

                    parents.put(childNode, currentNode);

                }

            }

        }

    }



    public Alignment alignBrute() {

        Queue<LatticeNodeHolder> open = new PriorityQueue<>();

        Map<LatticeNode, Integer> distance = new HashMap<>();

        Map<LatticeNode, LatticeNode> parents = new HashMap<>();



        LatticeNode sourceNode = getSourceNode();

        LatticeNode targetNode = getTargetNode();



        open.add(new LatticeNodeHolder(sourceNode, 0));

        distance.put(sourceNode, 0);

        parents.put(sourceNode, null);



        Set<LatticeNode> closed = new HashSet<>();



        while (true) {

            LatticeNode currentNode = open.remove().getNode();



            if (currentNode.equals(targetNode)) {

                return tracebackPath(parents, distance.get(currentNode));

            }



            if (closed.contains(currentNode)) {

                continue;

            }



            closed.add(currentNode);



            for (LatticeNode childNode : currentNode.getChildren()) {

                if (closed.contains(childNode)) {

                    continue;

                }



                int tentativeCost = distance.get(currentNode) +

                                    getWeight(currentNode, childNode);

                Integer currentCost = distance.get(childNode);



                if (currentCost == null || currentCost > tentativeCost) {

                    open.add(new LatticeNodeHolder(childNode, tentativeCost));

                    distance.put(childNode, tentativeCost);

                    parents.put(childNode, currentNode);

                }

            }

        }

    }



    /**

     * Create the source node for the sequence alignment task, i.e., a lattice 

     * node with all coordinates set to zero.

     * 

     * @return the source node.

     */

    LatticeNode getSourceNode() {

        int sourceCoordinates = new int[sequenceArray.length];

        return new LatticeNode(this, sourceCoordinates);

    }



    /**

     * Create the target node for the sequence alignment task, i.e., a lattice

     * node with all coordinates set to the respective sequence lengths.

     * 

     * @return the target node.

     */

    LatticeNode getTargetNode() {

        int targetCoordinates = new int[sequenceArray.length];



        for (int i = 0; i != sequenceArray.length; ++i) {

            targetCoordinates[i] = sequenceArray[i].length();

        }



        return new LatticeNode(this, targetCoordinates);

    }



    Integer getWeight(LatticeNode tail,

                      LatticeNode head, 

                      int dimension1,

                      int dimension2) {

        int tailCoordinates = tail.getCoordinates();

        int headCoordinates = head.getCoordinates();



        if (tailCoordinates[dimension1] == headCoordinates[dimension1]) {

            //System.out.println("1");

            return gapPenalty;

        } else if (tailCoordinates[dimension2] == headCoordinates[dimension2]) {

            //System.out.println("2");

            return gapPenalty;

        } else {

//            System.out.println("3");

            char character1 = sequenceArray[dimension1]

                    .charAt(tailCoordinates[dimension1]);



            char character2 = sequenceArray[dimension2]

                    .charAt(tailCoordinates[dimension2]);

            return costMatrix.getCost(character1, character2);

        }

    }



    Integer getWeight(LatticeNode tail, LatticeNode head) {

        // Extract the column represented by taking a single hop from 'tail' to

        // 'head':

        int tailCoordinates = tail.getCoordinates();

        int headCoordinates = head.getCoordinates();



        for (int i = 0; i < sequenceArray.length; ++i) {

            if (tailCoordinates[i] + 1 == headCoordinates[i]) {

                column[i] = sequenceArray[i].charAt(tailCoordinates[i]);

            } else {

                column[i] = AminoAcidAlphabet.GAP_CHARACTER;

            }

        }



        // Compute the hop cost as the sum of pairwise hops in any plane:

        int cost = 0;



        for (int i = 0; i < column.length; ++i) {

            char character1 = column[i];



            for (int j = i + 1; j < column.length; ++j) {

                char character2 = column[j];



                if (character1 != AminoAcidAlphabet.GAP_CHARACTER) {

                    if (character2 != AminoAcidAlphabet.GAP_CHARACTER) {

                        cost += costMatrix.getCost(character1, character2);

                    } else {

                        cost += gapPenalty;

                    }

                } else {

                    // character1 IS the gap character:

                    if (character2 != AminoAcidAlphabet.GAP_CHARACTER) {

                        cost += gapPenalty;

                    } else {

                        // Do nothing since we have a pair (gap, gap).

                    }

                }

            }

        }



        return cost;

    }



    String getSequenceArray() {

        return sequenceArray;

    }



    private Alignment tracebackPath(Map<LatticeNode, LatticeNode> parents,

                                    Integer cost) {

        List<LatticeNode> path = new ArrayList<>();

        LatticeNode node = getTargetNode();



        while (node != null) {

            path.add(node);

            node = parents.get(node);

        }



        Collections.<LatticeNode>reverse(path);



        String strings = new String[getSequenceArray().length];

        StringBuilder stringBuilders = new StringBuilder[strings.length];



        for (int i = 0; i < stringBuilders.length; ++i) {

            stringBuilders[i] = new StringBuilder();

        }



        for (int i = 1; i < path.size(); ++i) {

            LatticeNode tail = path.get(i - 1);

            LatticeNode head = path.get(i);



            int tailCoordinates = tail.getCoordinates();

            int headCoordinates = head.getCoordinates();



            for (int j = 0; j < tailCoordinates.length; ++j) {

                if (tailCoordinates[j] != headCoordinates[j]) {

                    stringBuilders[j].

                            append(sequenceArray[j].charAt(tailCoordinates[j]));

                } else {

                    stringBuilders[j].append(AminoAcidAlphabet.GAP_CHARACTER);

                }

            }

        }



        for (int i = 0; i < strings.length; ++i) {

            strings[i] = stringBuilders[i].toString();

        }



        return new Alignment(strings, cost);

    }



    private String checkIsValidGenomicSequence(String string) {

        Set<Character> characterSet = 

                AminoAcidAlphabet.getAminoAcidAlphabet()

                                 .getCharacterSet();



        for (char c : string.toCharArray()) {

            if (!characterSet.contains(c)) {

                throw new IllegalArgumentException("Unknown amino acid: " + c);

            }

        }



        return string;

    }



    private static final class LatticeNodeHolder

            implements Comparable<LatticeNodeHolder> {



        private final Integer fScore;

        private final LatticeNode node;



        LatticeNodeHolder(LatticeNode node, Integer fScore) {

            this.fScore = fScore;

            this.node = node;

        }



        LatticeNode getNode() {

            return node;

        }



        @Override

        public int compareTo(LatticeNodeHolder o) {

            return Integer.compare(fScore, o.fScore);

        }

    }

}

HeuristicFunctionComputer.java

package net.coderodde.bio.msa;



import java.awt.Point;

import java.util.HashSet;

import java.util.PriorityQueue;

import java.util.Queue;

import java.util.Set;



final class HeuristicFunctionComputer {



    HeuristicFunction computeHeuristicFunction(

            MultipleSequenceAlignmentInstance instance) {

        String sequenceArray = instance.getSequenceArray();

        HeuristicFunction heuristicFunction = new HeuristicFunction(instance);



        for (int dimension1 = 0; 

                dimension1 < sequenceArray.length; 

                dimension1++) {

            for (int dimension2 = dimension1 + 1; 

                    dimension2 < sequenceArray.length; 

                    dimension2++) {

                loadPartialHeuristicFunction(heuristicFunction,

                                             dimension1,

                                             dimension2,

                                             instance);

            }

        }



        return heuristicFunction;

    }



    // Basically, this method runs Dijkstra backwards from the target node until

    // the entire 2D-grid is explored.

    private void loadPartialHeuristicFunction(

            HeuristicFunction heuristicFunction,

            int dimension1,

            int dimension2,

            MultipleSequenceAlignmentInstance instance) {

        LatticeNode target = instance.getTargetNode();

        Queue<LatticeNodeHolder> open = new PriorityQueue<>();

        Set<LatticeNode> closed = new HashSet<>();

        open.add(new LatticeNodeHolder(target, 0));

        Point point = new Point();



        while (!open.isEmpty()) {

            LatticeNodeHolder currentNodeHolder = open.remove();

            LatticeNode currentNode = currentNodeHolder.getLatticeNode();



            if (closed.contains(currentNode)) {

                continue;

            }



            closed.add(currentNode);



            Integer currentNodeCost = currentNodeHolder.getCost();

            extractPoint(point, currentNode, dimension1, dimension2);

            heuristicFunction.put(dimension1,

                                  dimension2, 

                                  new Point(point),

                                  currentNodeCost);



            for (LatticeNode parent : currentNode.getParents()) {

                if (closed.contains(parent)) {

                    continue;

                }



                int tentativeCost =

                        currentNodeCost + instance.getWeight(parent, 

                                                             currentNode, 

                                                             dimension1, 

                                                             dimension2);

                extractPoint(point, parent, dimension1, dimension2);



                if (!heuristicFunction.containsPartial(dimension1,

                                                       dimension2, 

                                                       point)

                        || heuristicFunction.getPartial(dimension1,

                                                        dimension2,

                                                        point)

                        > tentativeCost) {

                    open.add(new LatticeNodeHolder(parent, tentativeCost));

                }

            }

        }

    }



    private static void extractPoint(Point point,

                                     LatticeNode latticeNode, 

                                     int i, 

                                     int j) {

        int coordinates = latticeNode.getCoordinates();

        point.x = coordinates[i];

        point.y = coordinates[j];

    }



    private static final class LatticeNodeHolder 

         implements Comparable<LatticeNodeHolder> {



        private final LatticeNode node;

        private final Integer cost;



        LatticeNodeHolder(LatticeNode node, Integer cost) {

            this.node = node;

            this.cost = cost;

        }



        LatticeNode getLatticeNode() {

            return node;

        }



        Integer getCost() {

            return cost;

        }



        @Override

        public int compareTo(LatticeNodeHolder o) {

            return Integer.compare(cost, o.cost);

        }

    }

}

LatticeNode.java

package net.coderodde.bio.msa;



import com.sun.corba.se.spi.orb.OperationFactory;

import java.util.Arrays;



final class LatticeNode {



    /**

     * The problem instance this lattice node contributes to.

     */

    private final MultipleSequenceAlignmentInstance instance;



    /**

     * The coordinates of this node in the entire lattice.

     */

    private final int coordinates;



    LatticeNode(MultipleSequenceAlignmentInstance instance, int coordinates) {

        this.instance = instance;

        this.coordinates = coordinates;

    }



    @Override

    public boolean equals(Object o) {

        return Arrays.equals(coordinates, ((LatticeNode) o).coordinates);

    }



    // Generated by NetBeans 8.1:

    @Override

    public int hashCode() {

        int hash = 7;

        hash = 41 * hash + Arrays.hashCode(this.coordinates);

        return hash;

    }



    @Override

    public String toString() {

        StringBuilder sb = new StringBuilder("[");

        String separator = "";



        for (int coordinate : coordinates) {

            sb.append(separator).append(coordinate);

            separator = ", ";

        }



        return sb.append("]").toString();

    }



    LatticeNode getChildren() {

        // Find out in how many dimension we can move forward:

        int dimensionsNotReached = 0;

        String sequenceArray = instance.getSequenceArray();

        boolean inclusionMap = new boolean[coordinates.length];



        for (int i = 0; i < sequenceArray.length; ++i) {

            if (coordinates[i] < sequenceArray[i].length()) {

                dimensionsNotReached++;

                // We can make a step forward in the direction of ith dimension:

                inclusionMap[i] = true;

            }

        }



        // Create the array of children:

        int numberOfChildren = pow2(dimensionsNotReached) - 1;

        LatticeNode children = new LatticeNode[numberOfChildren];

        loadChildren(children, inclusionMap);



        // Convert offsets to actual child coordinates:

        for (LatticeNode child : children) {

            child.addOffsets(coordinates);

        }



        return children;

    }



    LatticeNode getParents() {

        // Find out in how many dimensions we can move BACKward:

        int dimensionsNotReached = 0;

        String sequenceArray = instance.getSequenceArray();

        boolean inclusionMap = new boolean[coordinates.length];



        for (int i = 0; i < sequenceArray.length; ++i) {

            if (coordinates[i] > 0) {

                dimensionsNotReached++;

                // We can make a step backwards in the direction of ith 

                // dimension:

                inclusionMap[i] = true;

            }

        }



        // Create the array of parents:

        int numberOfParents = pow2(dimensionsNotReached) - 1;

        LatticeNode parents = new LatticeNode[numberOfParents];

        loadParents(parents, inclusionMap);



        // Convert the offsets to actual parent coordinates:

        for (LatticeNode parent : parents) {

            parent.addOffsets(coordinates);

        }



        return parents;

    }



    int getCoordinates() {

        return coordinates;

    }



    private void loadChildren(LatticeNode children, boolean inclusionMap) {

        int coords = new int[this.coordinates.length];



        for (int i = 0; i != children.length; ++i) {

            increment(coords, inclusionMap);

            children[i] = new LatticeNode(instance, coords.clone());

        }

    }



    private void loadParents(LatticeNode parents, boolean inclusionMap) {

        int coords = new int[this.coordinates.length];



        for (int i = 0; i != parents.length; ++i) {

            decrement(coords, inclusionMap);

            parents[i] = new LatticeNode(instance, coords.clone());

        }

    }



    private static void increment(int coordinates, boolean inclusionMap) {

        for (int i = 0; i < coordinates.length; ++i) {

            if (!inclusionMap[i]) {

                continue;

            }



            if (coordinates[i] == 0) {

                coordinates[i] = 1;

                return;

            }



            coordinates[i] = 0;

        }

    }



    private static void decrement(int coordinates, boolean inclusionMap) {

        for (int i = 0; i < coordinates.length; ++i) {

            if (!inclusionMap[i]) {

                continue;

            }



            if (coordinates[i] == 0) {

                coordinates[i] = -1;

                return;

            }



            coordinates[i] = 0;

        }

    }



    private void addOffsets(int offsets) {

        for (int i = 0; i < coordinates.length; ++i) {

            coordinates[i] += offsets[i];

        }

    }



    private static int pow2(int exponent) {

        int ret = 1;



        for (int e = 0; e < exponent; ++e) {

            ret *= 2;

        }



        return ret;

    }

}

PAM250CostMatrix.java

package net.coderodde.bio.msa;



import java.util.HashMap;

import java.util.Map;

import net.coderodde.bio.msa.AminoAcidAlphabet;

import net.coderodde.bio.msa.CostMatrix;



public final class PAM250CostMatrix implements CostMatrix<Integer> {



    private static PAM250CostMatrix instance;



    private final Map<Character, Map<Character, Integer>> m = new HashMap<>();



    public static PAM250CostMatrix getPAM250CostMatrix() {

        if (instance == null) {

            instance = new PAM250CostMatrix();

        }



        return instance;

    } 



    private PAM250CostMatrix() {

        AminoAcidAlphabet alphabet = AminoAcidAlphabet.getAminoAcidAlphabet();



        alphabet.getCharacterSet().stream().forEach((character) -> {

            m.put(character, new HashMap<>());

        });



        m.get('A').put('A', -2);

        m.get('R').put('A', 2);

        m.get('A').put('R', 2);

        m.get('R').put('R', -6);

        m.get('N').put('A', 0);

        m.get('A').put('N', 0);

        m.get('N').put('R', 0);

        m.get('R').put('N', 0);

        m.get('N').put('N', -2);

        m.get('D').put('A', 0);

        m.get('A').put('D', 0);

        m.get('D').put('R', 1);

        m.get('R').put('D', 1);

        m.get('D').put('N', -2);

        m.get('N').put('D', -2);

        m.get('D').put('D', -4);

        m.get('C').put('A', 2);

        m.get('A').put('C', 2);

        m.get('C').put('R', 4);

        m.get('R').put('C', 4);

        m.get('C').put('N', 4);

        m.get('N').put('C', 4);

        m.get('C').put('D', 5);

        m.get('D').put('C', 5);

        m.get('C').put('C', -4);

        m.get('Q').put('A', 0);

        m.get('A').put('Q', 0);

        m.get('Q').put('R', -1);

        m.get('R').put('Q', -1);

        m.get('Q').put('N', -1);

        m.get('N').put('Q', -1);

        m.get('Q').put('D', -2);

        m.get('D').put('Q', -2);

        m.get('Q').put('C', 5);

        m.get('C').put('Q', 5);

        m.get('Q').put('Q', -4);

        m.get('E').put('A', 0);

        m.get('A').put('E', 0);

        m.get('E').put('R', 1);

        m.get('R').put('E', 1);

        m.get('E').put('N', -1);

        m.get('N').put('E', -1);

        m.get('E').put('D', -3);

        m.get('D').put('E', -3);

        m.get('E').put('C', 5);

        m.get('C').put('E', 5);

        m.get('E').put('Q', -2);

        m.get('Q').put('E', -2);

        m.get('E').put('E', -4);

        m.get('G').put('A', -1);

        m.get('A').put('G', -1);

        m.get('G').put('R', 3);

        m.get('R').put('G', 3);

        m.get('G').put('N', 0);

        m.get('N').put('G', 0);

        m.get('G').put('D', -1);

        m.get('D').put('G', -1);

        m.get('G').put('C', 3);

        m.get('C').put('G', 3);

        m.get('G').put('Q', 1);

        m.get('Q').put('G', 1);

        m.get('G').put('E', 0);

        m.get('E').put('G', 0);

        m.get('G').put('G', -5);

        m.get('H').put('A', 1);

        m.get('A').put('H', 1);

        m.get('H').put('R', -2);

        m.get('R').put('H', -2);

        m.get('H').put('N', -2);

        m.get('N').put('H', -2);

        m.get('H').put('D', -1);

        m.get('D').put('H', -1);

        m.get('H').put('C', 3);

        m.get('C').put('H', 3);

        m.get('H').put('Q', -3);

        m.get('Q').put('H', -3);

        m.get('H').put('E', -1);

        m.get('E').put('H', -1);

        m.get('H').put('G', 2);

        m.get('G').put('H', 2);

        m.get('H').put('H', -6);

        m.get('I').put('A', 1);

        m.get('A').put('I', 1);

        m.get('I').put('R', 2);

        m.get('R').put('I', 2);

        m.get('I').put('N', 2);

        m.get('N').put('I', 2);

        m.get('I').put('D', 2);

        m.get('D').put('I', 2);

        m.get('I').put('C', 2);

        m.get('C').put('I', 2);

        m.get('I').put('Q', 2);

        m.get('Q').put('I', 2);

        m.get('I').put('E', 2);

        m.get('E').put('I', 2);

        m.get('I').put('G', 3);

        m.get('G').put('I', 3);

        m.get('I').put('H', 2);

        m.get('H').put('I', 2);

        m.get('I').put('I', -5);

        m.get('L').put('A', 2);

        m.get('A').put('L', 2);

        m.get('L').put('R', 3);

        m.get('R').put('L', 3);

        m.get('L').put('N', 3);

        m.get('N').put('L', 3);

        m.get('L').put('D', 4);

        m.get('D').put('L', 4);

        m.get('L').put('C', 6);

        m.get('C').put('L', 6);

        m.get('L').put('Q', 2);

        m.get('Q').put('L', 2);

        m.get('L').put('E', 3);

        m.get('E').put('L', 3);

        m.get('L').put('G', 4);

        m.get('G').put('L', 4);

        m.get('L').put('H', 2);

        m.get('H').put('L', 2);

        m.get('L').put('I', -2);

        m.get('I').put('L', -2);

        m.get('L').put('L', -6);

        m.get('K').put('A', 1);

        m.get('A').put('K', 1);

        m.get('K').put('R', -3);

        m.get('R').put('K', -3);

        m.get('K').put('N', -1);

        m.get('N').put('K', -1);

        m.get('K').put('D', 0);

        m.get('D').put('K', 0);

        m.get('K').put('C', 5);

        m.get('C').put('K', 5);

        m.get('K').put('Q', -1);

        m.get('Q').put('K', -1);

        m.get('K').put('E', 0);

        m.get('E').put('K', 0);

        m.get('K').put('G', 2);

        m.get('G').put('K', 2);

        m.get('K').put('H', 0);

        m.get('H').put('K', 0);

        m.get('K').put('I', 2);

        m.get('I').put('K', 2);

        m.get('K').put('L', 3);

        m.get('L').put('K', 3);

        m.get('K').put('K', -5);

        m.get('M').put('A', 1);

        m.get('A').put('M', 1);

        m.get('M').put('R', 0);

        m.get('R').put('M', 0);

        m.get('M').put('N', 2);

        m.get('N').put('M', 2);

        m.get('M').put('D', 3);

        m.get('D').put('M', 3);

        m.get('M').put('C', 5);

        m.get('C').put('M', 5);

        m.get('M').put('Q', 1);

        m.get('Q').put('M', 1);

        m.get('M').put('E', 2);

        m.get('E').put('M', 2);

        m.get('M').put('G', 3);

        m.get('G').put('M', 3);

        m.get('M').put('H', 2);

        m.get('H').put('M', 2);

        m.get('M').put('I', -2);

        m.get('I').put('M', -2);

        m.get('M').put('L', -4);

        m.get('L').put('M', -4);

        m.get('M').put('K', 0);

        m.get('K').put('M', 0);

        m.get('M').put('M', -6);

        m.get('F').put('A', 4);

        m.get('A').put('F', 4);

        m.get('F').put('R', 4);

        m.get('R').put('F', 4);

        m.get('F').put('N', 4);

        m.get('N').put('F', 4);

        m.get('F').put('D', 6);

        m.get('D').put('F', 6);

        m.get('F').put('C', 4);

        m.get('C').put('F', 4);

        m.get('F').put('Q', 5);

        m.get('Q').put('F', 5);

        m.get('F').put('E', 5);

        m.get('E').put('F', 5);

        m.get('F').put('G', 5);

        m.get('G').put('F', 5);

        m.get('F').put('H', 2);

        m.get('H').put('F', 2);

        m.get('F').put('I', -1);

        m.get('I').put('F', -1);

        m.get('F').put('L', -2);

        m.get('L').put('F', -2);

        m.get('F').put('K', 5);

        m.get('K').put('F', 5);

        m.get('F').put('M', 0);

        m.get('M').put('F', 0);

        m.get('F').put('F', -9);

        m.get('P').put('A', -1);

        m.get('A').put('P', -1);

        m.get('P').put('R', 0);

        m.get('R').put('P', 0);

        m.get('P').put('N', 1);

        m.get('N').put('P', 1);

        m.get('P').put('D', 1);

        m.get('D').put('P', 1);

        m.get('P').put('C', 3);

        m.get('C').put('P', 3);

        m.get('P').put('Q', 0);

        m.get('Q').put('P', 0);

        m.get('P').put('E', 1);

        m.get('E').put('P', 1);

        m.get('P').put('G', 1);

        m.get('G').put('P', 1);

        m.get('P').put('H', 0);

        m.get('H').put('P', 0);

        m.get('P').put('I', 2);

        m.get('I').put('P', 2);

        m.get('P').put('L', 3);

        m.get('L').put('P', 3);

        m.get('P').put('K', 1);

        m.get('K').put('P', 1);

        m.get('P').put('M', 2);

        m.get('M').put('P', 2);

        m.get('P').put('F', 5);

        m.get('F').put('P', 5);

        m.get('P').put('P', -6);

        m.get('S').put('A', -1);

        m.get('A').put('S', -1);

        m.get('S').put('R', 0);

        m.get('R').put('S', 0);

        m.get('S').put('N', -1);

        m.get('N').put('S', -1);

        m.get('S').put('D', 0);

        m.get('D').put('S', 0);

        m.get('S').put('C', 0);

        m.get('C').put('S', 0);

        m.get('S').put('Q', 1);

        m.get('Q').put('S', 1);

        m.get('S').put('E', 0);

        m.get('E').put('S', 0);

        m.get('S').put('G', -1);

        m.get('G').put('S', -1);

        m.get('S').put('H', 1);

        m.get('H').put('S', 1);

        m.get('S').put('I', 1);

        m.get('I').put('S', 1);

        m.get('S').put('L', 3);

        m.get('L').put('S', 3);

        m.get('S').put('K', 0);

        m.get('K').put('S', 0);

        m.get('S').put('M', 2);

        m.get('M').put('S', 2);

        m.get('S').put('F', 3);

        m.get('F').put('S', 3);

        m.get('S').put('P', -1);

        m.get('P').put('S', -1);

        m.get('S').put('S', -3);

        m.get('T').put('A', -1);

        m.get('A').put('T', -1);

        m.get('T').put('R', 1);

        m.get('R').put('T', 1);

        m.get('T').put('N', 0);

        m.get('N').put('T', 0);

        m.get('T').put('D', 0);

        m.get('D').put('T', 0);

        m.get('T').put('C', 2);

        m.get('C').put('T', 2);

        m.get('T').put('Q', 1);

        m.get('Q').put('T', 1);

        m.get('T').put('E', 0);

        m.get('E').put('T', 0);

        m.get('T').put('G', 0);

        m.get('G').put('T', 0);

        m.get('T').put('H', 1);

        m.get('H').put('T', 1);

        m.get('T').put('I', 0);

        m.get('I').put('T', 0);

        m.get('T').put('L', 2);

        m.get('L').put('T', 2);

        m.get('T').put('K', 0);

        m.get('K').put('T', 0);

        m.get('T').put('M', 1);

        m.get('M').put('T', 1);

        m.get('T').put('F', 2);

        m.get('F').put('T', 2);

        m.get('T').put('P', 0);

        m.get('P').put('T', 0);

        m.get('T').put('S', -1);

        m.get('S').put('T', -1);

        m.get('T').put('T', -3);

        m.get('W').put('A', 6);

        m.get('A').put('W', 6);

        m.get('W').put('R', -2);

        m.get('R').put('W', -2);

        m.get('W').put('N', 4);

        m.get('N').put('W', 4);

        m.get('W').put('D', 7);

        m.get('D').put('W', 7);

        m.get('W').put('C', 8);

        m.get('C').put('W', 8);

        m.get('W').put('Q', 5);

        m.get('Q').put('W', 5);

        m.get('W').put('E', 7);

        m.get('E').put('W', 7);

        m.get('W').put('G', 7);

        m.get('G').put('W', 7);

        m.get('W').put('H', 3);

        m.get('H').put('W', 3);

        m.get('W').put('I', 5);

        m.get('I').put('W', 5);

        m.get('W').put('L', 2);

        m.get('L').put('W', 2);

        m.get('W').put('K', 3);

        m.get('K').put('W', 3);

        m.get('W').put('M', 4);

        m.get('M').put('W', 4);

        m.get('W').put('F', 0);

        m.get('F').put('W', 0);

        m.get('W').put('P', 6);

        m.get('P').put('W', 6);

        m.get('W').put('S', 2);

        m.get('S').put('W', 2);

        m.get('W').put('T', 5);

        m.get('T').put('W', 5);

        m.get('W').put('W', -17);

        m.get('Y').put('A', 3);

        m.get('A').put('Y', 3);

        m.get('Y').put('R', 4);

        m.get('R').put('Y', 4);

        m.get('Y').put('N', 2);

        m.get('N').put('Y', 2);

        m.get('Y').put('D', 4);

        m.get('D').put('Y', 4);

        m.get('Y').put('C', 0);

        m.get('C').put('Y', 0);

        m.get('Y').put('Q', 4);

        m.get('Q').put('Y', 4);

        m.get('Y').put('E', 4);

        m.get('E').put('Y', 4);

        m.get('Y').put('G', 5);

        m.get('G').put('Y', 5);

        m.get('Y').put('H', 0);

        m.get('H').put('Y', 0);

        m.get('Y').put('I', 1);

        m.get('I').put('Y', 1);

        m.get('Y').put('L', 1);

        m.get('L').put('Y', 1);

        m.get('Y').put('K', 4);

        m.get('K').put('Y', 4);

        m.get('Y').put('M', 2);

        m.get('M').put('Y', 2);

        m.get('Y').put('F', -7);

        m.get('F').put('Y', -7);

        m.get('Y').put('P', 5);

        m.get('P').put('Y', 5);

        m.get('Y').put('S', 3);

        m.get('S').put('Y', 3);

        m.get('Y').put('T', 3);

        m.get('T').put('Y', 3);

        m.get('Y').put('W', 0);

        m.get('W').put('Y', 0);

        m.get('Y').put('Y', -10);

        m.get('V').put('A', 0);

        m.get('A').put('V', 0);

        m.get('V').put('R', 2);

        m.get('R').put('V', 2);

        m.get('V').put('N', 2);

        m.get('N').put('V', 2);

        m.get('V').put('D', 2);

        m.get('D').put('V', 2);

        m.get('V').put('C', 2);

        m.get('C').put('V', 2);

        m.get('V').put('Q', 2);

        m.get('Q').put('V', 2);

        m.get('V').put('E', 2);

        m.get('E').put('V', 2);

        m.get('V').put('G', 1);

        m.get('G').put('V', 1);

        m.get('V').put('H', 2);

        m.get('H').put('V', 2);

        m.get('V').put('I', -4);

        m.get('I').put('V', -4);

        m.get('V').put('L', -2);

        m.get('L').put('V', -2);

        m.get('V').put('K', 2);

        m.get('K').put('V', 2);

        m.get('V').put('M', -2);

        m.get('M').put('V', -2);

        m.get('V').put('F', 1);

        m.get('F').put('V', 1);

        m.get('V').put('P', 1);

        m.get('P').put('V', 1);

        m.get('V').put('S', 1);

        m.get('S').put('V', 1);

        m.get('V').put('T', 0);

        m.get('T').put('V', 0);

        m.get('V').put('W', 6);

        m.get('W').put('V', 6);

        m.get('V').put('Y', 2);

        m.get('Y').put('V', 2);

        m.get('V').put('V', -4);

    }



    @Override

    public Integer getCost(Character aminoAcidChar1, Character aminoAcidChar2) {

        try {

            return m.get(aminoAcidChar1).get(aminoAcidChar2);

        } catch (NullPointerException ex) {

            throw new IllegalArgumentException("Bad arguments: (" +

                    aminoAcidChar1 + ", " + aminoAcidChar2 + ")",

                    ex);

        }

    }

}

AminoAcidAlphabet.java

package net.coderodde.bio.msa;



import java.util.Collections;

import java.util.HashSet;

import java.util.Set;



public final class AminoAcidAlphabet {



    public static final Character GAP_CHARACTER = '-';

    private static AminoAcidAlphabet instance;



    private final Set<Character> alphabet = new HashSet<>();



    public static AminoAcidAlphabet getAminoAcidAlphabet() {

        if (instance == null) {

            instance = new AminoAcidAlphabet();

        }



        return instance;

    }



    private AminoAcidAlphabet() {

        String aminoAcids = "ACDEF" +

                            "GHIKL" +

                            "MNPQR" +

                            "STVWY";



        for (char c : aminoAcids.toCharArray()) {

            alphabet.add(c);

        }

    }



    public Set<Character> getCharacterSet() {

        return Collections.<Character>unmodifiableSet(alphabet);

    }

}

Alignment.java

package net.coderodde.bio.msa;



public class Alignment {



    private final String alignment;

    private final int cost;



    Alignment(String alignment, int cost) {

        this.alignment = alignment;

        this.cost = cost;

    }



    public String getAlignemnt() {

        return alignment;

    }



    public int getCost() {

        return cost;

    }



    @Override

    public String toString() {

        StringBuilder sb = new StringBuilder();

        String separator = "";



        for (String row : alignment) {

            sb.append(separator).append(row);

            separator = "n";

        }



        sb.append("nCost: ").append(cost);

        return sb.toString();

    }

}

CostMatrix.java

package net.coderodde.bio.msa;



public interface CostMatrix<C> {



    public C getCost(Character aminoAcidChar1, Character aminoAcidChar2);

}

HeuristicFunction.java

package net.coderodde.bio.msa;



import java.awt.Point;

import java.util.HashMap;

import java.util.Map;



final class HeuristicFunction {



    private final Map<Point, Integer> matrix;



    HeuristicFunction(MultipleSequenceAlignmentInstance instance) {

        int sequences = instance.getSequenceArray().length;

        this.matrix = new Map[sequences];



        for (int i = 0; i < sequences; ++i) {

            this.matrix[i] = new Map[sequences];



            for (int j = i + 1; j < sequences; ++j) {

                this.matrix[i][j] = new HashMap<>();

            }

        }

    }



    void put(int dimension1, int dimension2, Point point, Integer cost) {

        matrix[dimension1][dimension2].put(point, cost);

    }



    boolean containsPartial(int dimension1, int dimension2, Point point) {

        return matrix[dimension1][dimension2].containsKey(point);

    }



    Integer getPartial(int dimension1, int dimension2, Point point) {

        return matrix[dimension1][dimension2].get(point);

    }



    Integer get(LatticeNode node) {

        int cost = 0;

        int coordinates = node.getCoordinates();

        Point point = new Point();



        for (int dimension1 = 0; 

                dimension1 < coordinates.length; 

                dimension1++) {

            point.x = dimension1;



            for (int dimension2 = dimension1 + 1;

                    dimension2 < coordinates.length; 

                    dimension2++) {

                point.y = dimension2;

                cost += matrix[dimension1][dimension2].get(point);

            }

        }



        return cost;

    }

}

App.java

package net.coderodde.bio.msa;



final class App {



    private static final String SEQUENCES = {

        "ACGHKGMNPFQEKKFKLMNRW",

        "CFGPQWYRTLMEKKFKNR",

        "EACLMNRPQWTR",

        "TIMWAYHTMGIEKKFK"

    };



    public static void main(String args) {

        MultipleSequenceAlignmentInstance instance = 

                new MultipleSequenceAlignmentInstance(

                        PAM250CostMatrix.getPAM250CostMatrix(),

                        4, 

                        SEQUENCES);



        long start = System.currentTimeMillis();

        Alignment alignment1 = instance.align();

        long end = System.currentTimeMillis();



        System.out.println(alignment1);

        System.out.println(end - start + " ms.");

        System.out.println();



        start = System.currentTimeMillis();

        Alignment alignment2 = instance.alignBrute();

        end = System.currentTimeMillis();



        System.out.println(alignment2);

        System.out.println(end - start + " ms.");

    }

}

Output

The output for the demo is

Computed heuristic function in 1003 milliseconds.

-ACG-HKGMNPF-QEKKFKLMNRW

-CFGPQW-YRTL-MEKKFK--NR-

EAC--LM-NRP--Q-WT-R-----

-TI--MWAYHTMGIEKKFK-----

Cost: 53

1126 ms.



-ACG-HKGMNPF-QEKKFKLMNRW

-CFGPQW-YRTL-MEKKFK--NR-

EAC--LM-NRP--Q-WT-R-----

-TI--MWAYHTMGIEKKFK-----

Cost: 53

60 ms.

Critique request

As always, please tell me anything that comes to mind.

In AminoAcidAlphabet:

• 
  getCharacterSet always return same set, but you always create new instance.

In PAM250CostMatrix

• in constructor you can extract method to write add('A', 'A', -2) instead of m.get('A').put('A', -2)
  

In LatticeNode

• you can make hashCode simply do return Arrays.hashCode(this.coordinates);