# Multiple sequence alignment in Java

(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();

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;
}

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) {
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();

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;
}

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) {
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,
int dimension1,
int dimension2) {
int[] tailCoordinates = tail.getCoordinates();

//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
int[] tailCoordinates = tail.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) {
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);

int[] tailCoordinates = tail.getCoordinates();

for (int j = 0; j < tailCoordinates.length; ++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++) {
dimension1,
dimension2,
instance);
}
}

return heuristicFunction;
}

// Basically, this method runs Dijkstra backwards from the target node until
// the entire 2D-grid is explored.
HeuristicFunction heuristicFunction,
int dimension1,
int dimension2,
MultipleSequenceAlignmentInstance instance) {
LatticeNode target = instance.getTargetNode();
Queue<LatticeNodeHolder> open = new PriorityQueue<>();
Set<LatticeNode> closed = new HashSet<>();
Point point = new Point();

while (!open.isEmpty()) {
LatticeNodeHolder currentNodeHolder = open.remove();
LatticeNode currentNode = currentNodeHolder.getLatticeNode();

if (closed.contains(currentNode)) {
continue;
}

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) {
}
}
}
}

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];

// Convert offsets to actual child coordinates:
for (LatticeNode child : children) {
}

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];

// Convert the offsets to actual parent coordinates:
for (LatticeNode parent : parents) {
}

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;
}
}

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()) {
}
}

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);