Optimizing Needleman-Wunsch algorithm in Java

Question

For a class assignment, I was required to implement the Needleman-Wunsch algorithm in Java. I have already completed the assignment but I was wondering if there were any algorithmic (or even syntactic) improvements that could be made over my existing code to make it more efficient.

public class DNASequence {
    protected final String seq_1, seq_2;    // The two sequences to be analyzed
    public int alignmentScore;              // Using Needleman-Wunsch algorithm
    protected Node[][] matrix;              // Store scores and indels

    protected final int matchScore, mismatchScore, indel;

    // Strings to be used to print the DNA sequnce analysis
    String top = "";        // Sequence 1
    String buffer = "";     // Matches between Sequnce 1 & 2
    String bottom = "";     // Sequence 2

    public DNASequence(String s1, String s2, ScoreScheme s) {
        // I use a ■ as a buffer so that the Sequence string aligns properly
        // with the indels and scores within the matrix.
        seq_1 = "\u25A0" + s1;
        seq_2 = "\u25A0" + s2;

        // Setup the scoring schema. For this program, I am implementing a simple
        // gap cost instead of complicating it further with gap extention and
        // opening costs.
        matchScore = s.matchScore;
        mismatchScore = s.mismatchScore;
        indel = s.indel;

        // I instiate the matrix and progressively build the indels on the first
        // row and the first column.
        matrix = new Node[seq_1.length()][seq_2.length()];
        for (int i = 0; i < seq_1.length(); i++)
            matrix[i][0] = new Node(i * indel, i, 0);

        for (int i = 0; i < seq_2.length(); i++)
            matrix[0][i] = new Node(i * indel, 0, i);
    }

    // Helper method that helps decide what kind of match/mismatch score to use.
    protected int score(int i, int j) {
        if (seq_1.charAt(i) == seq_2.charAt(j))
            return matchScore;
        else
            return mismatchScore;
    }

    // Helper method that implements the Needleman-Wunsch algo on a local level.
    protected Node match(int i, int j) {
        Node s1,s2,s3;
        s1 = new Node(matrix[i-1][j-1].score + score(i, j), i, j);
        s2 = new Node(matrix[i-1][j].score + indel, i, j);
        s3 = new Node(matrix[i][j-1].score + indel, i, j);

        // Since the aim of the Needleman-Wunsch algo is to find the shortest path
        // with the lowest cost and highest score, I check whichever of the viable
        // nodes nets me the highest score (up till that point). Once found, I set
        // a pointer back to the previous node (to make it easier to traceback)
        // and return the current node (to be placed in the matrix).
        Node largest = new Node(Math.max(s1.score, Math.max(s2.score, s3.score)), i, j);
        if (s1.compareTo(largest) == 0)
            largest.prev = matrix[i-1][j-1];
        else if(s2.compareTo(largest) == 0)
            largest.prev = matrix[i-1][j];
        else
            largest.prev = matrix[i][j-1];

        return largest;
    }

    // Runs the Needleman-Wunsch algo on every node in the matrix, sets the aligment
    // score and returns the last node in the matrix -- the one that holds the score.
    public Node runAnalysis() {
        for (int i = 1; i < seq_1.length(); i++) {
            for (int j = 1; j < seq_2.length(); j++){
                matrix[i][j] = match(i, j);
            }
        }
        alignmentScore = matrix[seq_1.length()-1][seq_2.length()-1].score;
        return matrix[seq_1.length()-1][seq_2.length()-1];
    }

    // Helper method that progressively builds the analysis result. It returns the 
    // 'tail' because we may still need to do some work on it.
    protected Node traceHelper(Node curr) {
        while (curr.prev != null) {
            // Print out the 'path'
            // System.out.print(curr.score + "[" + curr.i + ", " + curr.j + "] -> ");

            if (curr.i - curr.prev.i == 1 && curr.j - curr.prev.j == 1){    // If the path leads diagonal
                boolean x = seq_1.charAt(curr.i) == seq_2.charAt(curr.j) ? true : false;
                // System.out.println("Going diag: " + x);
                if(x){
                    top = seq_1.charAt(curr.i) + " " + top;
                    buffer = "|" + " " + buffer;
                    bottom = seq_2.charAt(curr.j) + " " + bottom;
                }else{
                    top = seq_1.charAt(curr.i) + " " + top;
                    buffer = " " + " " + buffer;
                    bottom = seq_2.charAt(curr.j) + " " + bottom;
                }
            }else if (curr.i - curr.prev.i == 1){                           // If the path leads up
                // System.out.println("Going up: " + indel);
                top = seq_1.charAt(curr.i) + " " + top;
                buffer = " " + " " + buffer;
                bottom = "-" + " " + bottom;                                // If the path leads left
            }else if (curr.j - curr.prev.j == 1){
                // System.out.println("Going left: " + indel);
                top = "-" + " " + top;
                buffer = " " + " " + buffer;
                bottom = seq_2.charAt(curr.j) + " " + bottom;
            }

            curr = curr.prev;
        }

        return curr;
    }

    // Traceback from the last node in the matrix back to the first indel node.
    public void traceback() {
        Node curr = matrix[seq_1.length()-1][seq_2.length()-1];

        curr = traceHelper(curr);

        // Sometimes the tail or the traceback path ends on an indel. As a result,
        // we end up with an incomplete path because the algorithm doesn't see past
        // it. To counter that, I check if I am on an indel and if so, I help it
        // move along the path back to 0,0.
        while (curr.i != 0 || curr.j != 0) {
            if (curr.i != 0 && curr.j == 0){
                curr.prev = matrix[curr.i-1][curr.j];
                curr = traceHelper(curr);
            }else if (curr.i == 0 && curr.j != 0) {
                curr.prev = matrix[curr.i][curr.j-1];
                curr = traceHelper(curr);
            }
        }

        // Print out the DNA sequence analysis
        System.out.println(top);
        System.out.println(buffer);
        System.out.println(bottom);
    }

    // Print the matrix.
    public void printMatrix() {
        System.out.printf("%4s", "\u25A0");
        for (int i = 0; i < matrix[0].length; i++) {
            System.out.printf("%4s", seq_2.charAt(i));
        }
        System.out.println();
        for (int i = 0; i < matrix.length; i++) {
            System.out.printf("%4s", seq_1.charAt(i));
            for (int j = 0; j < matrix[i].length; j++) {
                System.out.printf("%4d", matrix[i][j].score);
            }
            System.out.println();
        }
    }

    // Creates a random DNA sequence -- for testing purposes.
    public static String randseq(int n) {
        char[] S = new char[n];
        String DNA = "ACGT";

        for (int i = 0; i < n; i++) {
            int r = (int)(Math.random() * 4);
            S[i] = DNA.charAt(r);
        }

        return new String(S);
    }

    public static void main(String[] args) {
        // Test with randomly generated DNA sequences
        String seq_1 = randseq(34);
        String seq_2 = randseq(32);

        // Or test with the following
        // String seq_1 = "CATTAATTACACTCTCGCACTCACCACCAAACATCCTAAACCCAGACAGGCCTCGACTCC";
        // String seq_2 = "ACTAAACAAGACTCGCCTGTCTAACTAGGGAGTTTATAATGAACCGTGGCGTAGACCA";

        ScoreScheme s = new ScoreScheme(2, -1, -2);             // MatchScore: 2; MismatchScore: -2; Indel: -2
        DNASequence dna = new DNASequence(seq_1, seq_2, s);     // Create a new DNA Sequence using two strands

        dna.runAnalysis();

        dna.traceback();
        System.out.println("The alignment score is " + dna.alignmentScore);

        // dna.printMatrix();
    }
}

class Node implements Comparable<Node>{
    int i, j;
    int score;
    Node prev;

    public Node(int score, int x, int y) {
        this.i = x;
        this.j = y;
        this.score = score;
        this.prev = null;
    }

    public int compareTo(Node n) {
        return this.score - n.score;
    }

    public String toString() {
        return ""+score;
    }
}

class ScoreScheme {
    int matchScore, mismatchScore, indel;

    public ScoreScheme(int m1, int m2, int i) {
        matchScore = m1;
        mismatchScore = m2;
        indel = i;
    }
}

Don't assume that everyone knows what you are talking about. What is the "Needleman-Wunsch algorithm"? — Simon Forsberg
– Simon Forsberg, Commented Dec 12, 2017 at 20:20
My apologies. The Needleman-Wunsch algorithm is used mostly in the field of bioinformatics to align protein sequences. It is also a great example of dynamic programming, hence its implementation in a CS class. Knowledge of the algorithm isn't necessarily required to understand the code hence why I left it out. — Najm Sheikh
– Najm Sheikh, Commented Dec 13, 2017 at 22:25

forsvarir · Accepted Answer · 2020-02-15 12:03:36Z

I'm not really going to address your core algorithm, however there are things that you can do to make your code more approachable.

Access modifiers

You've got an interesting collection of access modifiers used in your code. The only one that's missing is private. In general, you want the methods that will be interacted with from outside to be declared as public and everything else in your class to be declared as private. You should only be using protected if you're expecting there to be a derived implementation of the class which needs access to those members (typically this would be methods, not fields).

Constants can help

There's a few magic strings in your code that could really benefit from being constants. Not only does this offer the opportunity to give a descriptive name to the value, but it reduces that chances of you creating a typo/mismatch. So, for example:

private static final String MATRIX_PADDING_CHARACTER = "\u25A0";

Parameter names are important

Your ScoreScheme constructor looks like this:

public ScoreScheme(int m1, int m2, int i)

At the time you wrote the code, this probably made sense, without looking, does it still make sense? What's the difference between m1 and m2? It's easy enough to check by looking at the implementation, however the less context switching required, the easier it is to work with the code. Modern IDE's give helpful information, such as parameter names when you start coding a call to a method. If the parameters have sensible names, matchScore (m1) and mismatchScore (m2) are unlikely to get reversed by accident.

One declaration per line

Whilst the compiler allows you to declare multiple variables of the same type on the same line, I've rarely seen it done in practice.

int matchScore, mismatchScore, indel;

Having a separate line for each of the declarations doesn't cost a lot and can make the list easier to read. There's a very similar rule around single line if and for. Whilst you can do them, most places adopt an always include braces approach, because it's just too easy to make a mistake otherwise.

Scheme vs Members

You pass the ScoreScheme into your DNASequence class, however you then extract the values and store them in individual member variables.

matchScore = s.matchScore;
mismatchScore = s.mismatchScore;
indel = s.indel;

It's unclear why you wouldn't just store the ScoreScheme itself. This maintains the encapsulation of those variables so that it's clear they belong together.

Consistency

Again, this comes down to naming. The ScoreScheme members have a Score post-fix, apart from indel. They all represent scoring, either make them all '..Score' or not.

There's a similar mismatch in Node

this.i = x;
this.j = y;

The parameter is 'x', the field is 'i'. Picking a coordinate system and using it consistently throughout the code will make it easier to follow.

Dead Code

Dead / unused code should be removed before checkin / review. You've got some code that doesn't appear to ever be called (Node.toString) and other seconds of code that are just commented out, such as printlns that you used for debugging. The less noise there is in the code, the easier it is to focus on what it actually does.

Unnecessary checks

There's a redundant check in the else if, curr.j must not equal 0, or the while loop would already have terminated.

while (curr.i != 0 || curr.j != 0) {
  if (curr.i != 0 && curr.j == 0){
    curr.prev = matrix[curr.i-1][curr.j];
    curr = traceHelper(curr);
  }else if (curr.i == 0 && curr.j != 0) {
    curr.prev = matrix[curr.i][curr.j-1];
    curr = traceHelper(curr);
  }
}

As a side note, consistent spacing is important. Modern IDE's will auto-format your code so that it has consistent spacing. This makes it more approachable, but also helps to avoid false-positives when checking for differences in source control. If you always use the same format, then you don't need to worry about opening a file with a hundred changes in it only to find out all but one of them is somebody fixing the 'space before brace' style to their preference.

String vs StringBuilder

You have three buffers (top, buffer, bottom) which are all declared as strings. You build them up, starting from the right as you loop through the matrix. Typically, you'd be better off using a StringBuilder to perform these kinds of joining operations in a loop.

traceHelper and duplicate code

traceHelper is quite a big function, that has quite a lot of duplication in it. Essentially, you're prepending each new character and its match to the three buffers, with four scenarios. If I create a method:

final String NON_MATCHING_PADDING = " " + " ";

void insertPair(char seq1, boolean match, char seq2) {
    top.insert(0, seq1 + " ");
    buffer.insert(0, match ? "|" + " " : NON_MATCHING_PADDING );
    bottom.insert(0, seq2 + " ");
}

Then it can be simplified to:

protected Node traceHelper(Node curr) {
    while (curr.prev != null) {
        if (curr.i - curr.prev.i == 1 && curr.j - curr.prev.j == 1) {    // If the path leads diagonal
            insertPair(seq_1.charAt(curr.i), 
                       seq_1.charAt(curr.i) == seq_2.charAt(curr.j), 
                       seq_2.charAt(curr.j));
        } else if (curr.i - curr.prev.i == 1) {                           // If the path leads up
            insertPair(seq_1.charAt(curr.i), false, '-');
        } else if (curr.j - curr.prev.j == 1) {                           // If the path leads left
            insertPair('-', false, seq_2.charAt(curr.j));
        }
        curr = curr.prev;
    }

    return curr;
}

This makes the method a bit more approachable and highlights additional questions, such as could else if (curr.j - curr.prev.j == 1) be replaced by else, or is there another condition that needs to be handled?

Random

You only use if for generating test sequences, however prefer Random.nextInt over Math.random(). See discussion here.

Overall design

For me, it feels like there's a bit too much going on in the DNASequence class. It's responsible for taking to sequences and a score scheme, performing analysis, storing the results, outputting them to the console. Having an analyser that took in two sequences and a score scheme and produced some analysis results, which could then be printed would seem like a better split.

You also want to try to design your classes so that the caller needs to know as little as possible about how they work under the surface. If you've got two methods that are intrinsically linked, and need to be called in a specific order, then they should probably be a single public method on your class (you can have private methods if the split makes sense). At the moment, you have to call runAnalysis, followed by traceback to get the results. If you call traceback first, the application crashes. If you call it twice, you get double the result (because it doesn't clear the buffers before starting to append to them). From a client perspective, having them as separate methods doesn't really make sense, the traceback is really part of the analysis.

Stack Exchange Network

Optimizing Needleman-Wunsch algorithm in Java

1 Answer 1

You must log in to answer this question.

Hot Network Questions

Optimizing Needleman-Wunsch algorithm in Java

1 Answer 1

You must log in to answer this question.

Related

Hot Network Questions