动态编程之序列比对：Needleman-Wunsch 算法和Smith-Waterman算法

最新推荐文章于 2024-08-17 21:30:16 发布

HUGEHEADHUGE

最新推荐文章于 2024-08-17 21:30:16 发布

阅读量3.3w

点赞数

分类专栏： bioinformatics 文章标签：算法编程 string alignment matrix algorithm

bioinformatics 专栏收录该内容

4 篇文章 2 订阅

订阅专栏

动态编程之序列比对：Needleman-Wunsch 算法和Smith-Waterman算法

在比对两个序列时，不仅要考虑完全匹配的字符，还要考虑一个序列中的空格或间隙（或者，相反地，要考虑另一个序列中的插入部分）和不匹配，在序列比对中，需要找到最优的比对（最优比对大致是指要将匹配的数量最大化，将空格和不匹配的数量最小化）。如果要更正式些，您可以确定一个分数，为匹配的字符添加分数、为空格和不匹配的字符减去分数。

3. Algorithm for sequence alignment: dynamicprogramming

Making an alignment by hand is possible, but tedious. In some cases,when one has a lot of information about the proteins, such as activesite residues, secondary structure, 3D structure, mutations, etc, itmay still be necessary to make a manual alignment (or at least edit analignment) to fit all the data. The available automatic methods maynot be able to produce a good enough alignment in such cases.

Of course, we would like to have a completely automatic method toperform sequence alignment. The method of choice is based on so-calleddynamic programming, which is a general algorithm forsolving certain optimization problems. The word "programming" does notmean that it has to be a computer program; this is just mathematicaljargon for using a fixed set of rules to arrive at a solution.

For any automatic method to work, we need to be explicit about theassumptions that should go into it. We therefore need to have anexplicit scheme for the gap penalties and for the substitution matrix(see the previous page). The chosen gappenalties and substitution matrices are often collectively called thescoring scheme.

There are clearly many different possible scoring schemes. One mayalso complicate things further by allowing position-specific scores:If one knows from other sources (3D structure) that a gap should absolutely not be allowed in a certain part of a sequence, then thegap-open penalty could be set to a very high value in that part.

Given a scoring scheme, how does an alignment algorithm work? Let ususe the classical Needleman-Wunsch-Sellers algorithmto demonstrate how a dynamic-programming algorithm can work. Pleasenote that there are other variants of dynamic programming in sequenceanalysis.

The Needleman-Wunsch-Sellers algorithm sets up a matrix where eachsequence is placed along the sides of the matrix. Each element in thematrix represents the two residues of the sequences being aligned atthat position. To calculate the score in every position(i, j) one looks at the alignmentthat has already been made up to that point, and finds the best way to continue. Having gone through the entire matrix in this way, one cango back and trace which way through the matrix gives the best alignment.

Let us use the following gap-penalty function, wherek is the length of the gap,c_open the gap-open penalty constant, and c_length the gap-length penalty constant:

W(k) = c_open + c_length * k

The formulat describing the Needleman-Wunsch-Sellers method isrecursive, and for the position (i,j) is as follows, where D is valueof element (i, j) in the matrix and subst is the substitution matrix:

D_{i, j} = max {	D_{i - 1, j - 1} + subst(A_i, B_j)
	D_{i - 1, j - k} + W(k) (where k = 1, ..., j - 1)
	D_{i - k, j - 1} + W(k) (where k = 1, ..., i - 1)

After one has applied this to the matrix, one finds the optimalalignment by tracing backwards from the diagonal element backwards tothe previous highest value, and so on.

I have found a good tutorial describing dynamic programming forsequence alignment of the Needleman-Wunsch variant. The tutorialwas written by Eric C. Rouchka (Washington Universityin St. Louis), and walks through an example in detail.

Dynamic programming, a simplevariant, with no gap penalties and a simple substitution scoringscheme.
Dynamic programming, a moreadvanced variant, with gap penalties and a slightly morecomplicated substitution scoring scheme.

全局和局部序列比对:

全局序列比对尝试找到两个完整的序列 S1 和 S2 之间的最佳比对。如S1=GCCCTAGCG S2=GCGCAATG 如果设定每个匹配字符为1分，每个空格为-2分，每个不匹配为-1分，则下面的比对就是全局最优比对：S1'=GCCCTAGCG S2'=GCGC_AATG,连字符“_”代表空格。在 S2' 中有五个匹配字符，一个空格（或者反过来说，在 S1' 中有一个插入项），有三个不匹配字符。这样得到的分数是 (5×1) + (1×-2) + (3×-1) = 0，这是能够实现的最佳结果。

局部序列比对不必对两个完整的序列进行比对；可以在每个序列中使用某些部分来获得最大得分。使用同样的序列 S1 和 S2，以及同样的得分方案，可以得到以下局部最优比对 S1'' 和 S2''：S1=GCCCTAGCG S2=GCGCAATG S1''=GCG S2''=GCG，虽然这个局部比对恰好没有不匹配字符或空格，但是一般情况下，局部比对可能存在不匹配字符或空格。这个局部比对的得分是 (3×1) + (0×-2) + (0×-1) = 3。（最佳局部比对的得分要大于或等于最佳全局比对的得分，这是因为全局比对也属于局部比对)

Needleman-Wunsch 算法：

该算法用来计算全局比对。它的思路与 LCS 算法相似。这个算法也使用二维表格，一个序列沿顶部展开，一个序列沿左侧展开。而且也能通过以下三个途径到达每个单元格：1.来自上面的单元格，代表将左侧的字符与空格比对。2.来自左侧的单元格，代表将上面的字符与空格比对。3.来自左上侧的单元格，代表与左侧和上面的字符比对（可能匹配也可能不匹配）。该单元格的值来自于一下3个中的最大值：（1）上方的值-2 （2）左边的值-2 （3）如果该单元格所在的行于所在的列对应的字符相等，则为左上值加1，否则为左上值-1。

首先，必须初始化表格。这意味着填充第二行和第二列的分数和指针。填充第二行的操作意味着使用位于顶部的第一个序列中的字符，并使用空格，而不是使用左侧从上到下的序列中的第一个字符。空格的扣分是 -2，所以每次使用空格的时候，就给以前的单元格加了 -2 分。以前的单元格是左侧的单元格。这就说明了在第二行中为什么得到了 0, -2, -4, -6, ... 这样的序列。用相似的方式得到第二列的得分和指针.

接下来，需要填充余下的单元格。同 LCS 算法一样，对于每个单元格，都有三个选择，要从中选择最大的。可以从上面、左侧、左上侧到达每个单元格。假设 S1 和 S2 是要比对的字符串，S1' 和 S2' 是生成的比对中的字符串。从上面到达单元格相当于将左面的字符从 S2 加入 S2'，跳过上面的 S1 中的当前字符，并在 S1' 中加入一个空格。因为一个空格的分数是 -2，所以当前单元格的得分要从上面的单元格得分减 2 得到。类似的，将左边的单元格得分减 2，可以从左侧到达空单元格。

最后，可以将上面的字符加入到 S1' 中，将左边的字符加入到 S2' 中。这就相当于从左上侧进入空白单元格。这两个字符将会匹配，在这种情况下，新的得分就是左上侧单元格的得分减 1。在这三种可能性当中，选择得分最大的一个（如果得分相等，可以从得分高的单元格中从任选一个）。接下来，需要得到实际的比对字符串S1' 和S2'以及比对的得分。右下角单元格中的得分包含 S1 和 S2 的最大比对得分，就像在 LCS 算法中包含 LCS 的长度一样。而且，与 LCS 算法类似，要获得 S1' 和 S2'，要从右下角单元格开始沿着指针回溯，反向构建 S1' 和 S2'。从表格的构建过程可知，从上向下对应着将左侧字符从 S2 加入到 S2' 中，将空格加入 S1' 中；从左向右对应着将上面的字符从 S1 加入到 S1' 中，将空格加入 S2' 中；而向下和向右移动意味着分别将来自 S1 和 S2 的字符加入 S1' 和 S2'。

Smith-Waterman 算法：在 Smith-Waterman 算法中，不必比对整个序列。两个零长字符串即为得分为 0 的局部比对，这一事实表明在构建局部比对时，不需要使用负分。这样会造成进一步比对所得到的分数低于通过 “重设” 两个零长字符串所能得到的分数。而且，局部比对不需要到达任何一个序列的末端，所以也不需要从右下角开始回溯：可以从得分最高的单元格开始回溯。

这导致 Smith-Waterman 算法与 Needleman-Wunsch 算法存在着三个区别。

首先，在初始化阶段，第一行和第一列全填充为 0（而且第一行和第一列的指针均为空）。

第二，在填充表格时，如果某个得分为负，那么就用 0 代替，只对得分为正的单元格添加返回指针。

最后，在回溯的时候，从得分最高的单元格开始，回溯到得分为 0 的单元格为止。除此之外，回溯的方式与 Needleman-Wunsch 算法完全相同。

每个单元格的值来自于一下的最大值：：（1）上方的值-2 （2）左边的值-2 （3）如果该单元格所在的行于所在的列对应的字符相等，则为左上值加1，否则为左上值-1。但是切记，如果最大值为负数，则该单元格的值用0代替，不指向任何单元格。

两个算法的Java代码实现：

public abstract class SequenceAlignment extends DynamicProgramming {
    protected int match; //当匹配时的分值
    protected int mismatch; // 当不匹配时的分值
    protected int space; // 当为空格时的分值
    protected String[] alignments;
    public SequenceAlignment(String sequence1, String sequence2) {
        this(sequence1, sequence2, 1, -1, -1);
    }
    public SequenceAlignment(String sequence1, String sequence2, int match,
            int mismatch, int gap) {
        super(sequence1, sequence2);
        this.match = match;
        this.mismatch = mismatch;
        this.space = gap;
    }
    /*
     * (non-Javadoc) 两个算法通用的回溯
     *
     * @see com.ibm.combio.DynamicProgramming#getTracebace()
     */
    protected Object getTracebace() {
        StringBuffer align1Buf = new StringBuffer();
        StringBuffer align2Buf = new StringBuffer();
        Cell currentCell = getTracebackStartingCell();
        while (traceBackIsNotDone(currentCell)) {
            if (currentCell.getRow() - currentCell.getPrevCell().getRow() == 1) {
                align2Buf.insert(0, sequence2.charAt(currentCell.getRow() - 1));
            } else {
                align2Buf.insert(0, '-');
            }
            if (currentCell.getCol() - currentCell.getPrevCell().getCol() == 1) {
                align1Buf.insert(0, sequence1.charAt(currentCell.getCol() - 1));
            } else {
                align1Buf.insert(0, '-');
            }
            currentCell = currentCell.getPrevCell();
        }
        String[] alignments = new String[] { align1Buf.toString(),
                align2Buf.toString() };
        return alignments;
    }
    public int getAlignmentScore() {
        if (alignments == null) {
            getAlignment();
        }
        int score = 0;
        for (int i = 0; i < alignments[0].length(); i++) {
            char c1 = alignments[0].charAt(i);
            char c2 = alignments[1].charAt(i);
            if (c1 == '-' || c2 == '-') {
                score += space;
            } else if (c1 == c2) {
                score += match;
            } else {
                score += mismatch;
            }
        }
        return score;
    }
    public String[] getAlignment() {
        ensureTableIsFillIn();
        alignments = (String[]) getTracebace();
        return alignments;
    }
    protected abstract boolean traceBackIsNotDone(Cell currentCell); // 回溯停止判断
    protected abstract Cell getTracebackStartingCell(); // 取得回溯起始单元
}
2.NeedlemanWunsch 算法类
public class NeedlemanWunsch extends SequenceAlignment {
    public NeedlemanWunsch(String sequence1, String sequence2) {
        super(sequence1, sequence2);
    }
    public NeedlemanWunsch(String sequence1, String sequence2, int match,
            int mismatch, int gap) {
        super(sequence1, sequence2, match, mismatch, gap);
    }
    protected void fillInCell(Cell currentCell, Cell cellAbove,
            Cell cellToLeft, Cell cellAboveLeft) {
        int rowSpaceScore = cellAbove.getScore() + space;
        int colSpaceScore = cellToLeft.getScore() + space;
        int matchOrMismatchScore = cellAboveLeft.getScore();
        if (sequence2.charAt(currentCell.getRow() - 1) == sequence1
                .charAt(currentCell.getCol() - 1)) {
            matchOrMismatchScore += match;
        } else {
            matchOrMismatchScore += mismatch;
        }
        if (rowSpaceScore >= colSpaceScore) {
            if (matchOrMismatchScore >= rowSpaceScore) {
                currentCell.setScore(matchOrMismatchScore);
                currentCell.setPrevCell(cellAboveLeft);
            } else {
                currentCell.setScore(rowSpaceScore);
                currentCell.setPrevCell(cellAbove);
            }
        } else {
            if (matchOrMismatchScore >= colSpaceScore) {
                currentCell.setScore(matchOrMismatchScore);
                currentCell.setPrevCell(cellAboveLeft);
            } else {
                currentCell.setScore(colSpaceScore);
                currentCell.setPrevCell(cellToLeft);
            }
        }
    }
    @Override
    protected boolean traceBackIsNotDone(Cell currentCell) {
        return currentCell.getPrevCell() != null;
    }
    @Override
    protected Cell getTracebackStartingCell() {
        return scoreTable[scoreTable.length - 1][scoreTable[0].length - 1];
    }
    /*
     * (non-Javadoc)
     *
     * @see java.lang.Object#toString()
     */
    @Override
    public String toString() {
        return "[NeedlemanWunsch: sequence1=" + sequence1 + ", sequence2="
                + sequence2 + "]";
    }
    /* (non-Javadoc)
     * 初始化单元格的前驱指针
     * @see com.ibm.combio.DynamicProgramming#getInitialPointer(int, int)
     */
    protected Cell getInitialPointer(int row, int col) {
        if (row == 0 && col != 0) {
            return scoreTable[row][col - 1];
        } else if (col == 0 && row != 0) {
            return scoreTable[row - 1][col];
        } else {
            return null;
        }
    }
    /* (non-Javadoc)
     * 初始化单元格的值
     * @see com.ibm.combio.DynamicProgramming#getInitialScore(int, int)
     */
    protected int getInitialScore(int row, int col) {
        if (row == 0 && col != 0) {
            return col * space;
        } else if (col == 0 && row != 0) {
            return row * space;
        } else {
            return 0;
        }
    }
    @Override
    public Object getTracebace() {
        return super.getTracebace();
    }
}
3.SmithWaterman 类
public class SmithWaterman extends SequenceAlignment {
    private Cell highScoreCell; // 最高分值得单元
    public SmithWaterman(String sequence1, String sequence2) {
        super(sequence1, sequence2);
    }
    public SmithWaterman(String sequence1, String sequence2, int match,
            int mismatch, int gap) {
        super(sequence1, sequence2, match, mismatch, gap);
    }
    @Override
    protected void initialize() {
        // TODO Auto-generated method stub
        super.initialize();
        highScoreCell = scoreTable[0][0];
    }
    @Override
    protected Cell getTracebackStartingCell() {
        return highScoreCell;
    }
    @Override
    protected boolean traceBackIsNotDone(Cell currentCell) {
        return currentCell.getScore() != 0;
    }
    @Override
    public Object getTracebace() {
        return super.getTracebace();
    }
    @Override
    protected void fillInCell(Cell currentCell, Cell cellAbove,
            Cell cellToLeft, Cell cellLeftAbove) {
        int rowSpaceScore = cellAbove.getScore() + space;
        int colSpaceScore = cellToLeft.getScore() + space;
        int matchOrMismatchScore = cellLeftAbove.getScore();
        if (sequence1.charAt(currentCell.getCol() - 1) == sequence2
                .charAt(currentCell.getRow() - 1)) {
            matchOrMismatchScore += match;
        } else {
            matchOrMismatchScore += mismatch;
        }
        if (rowSpaceScore >= colSpaceScore) {
            if (matchOrMismatchScore >= rowSpaceScore) {
                if (matchOrMismatchScore > 0) {
                    currentCell.setScore(matchOrMismatchScore);
                    currentCell.setPrevCell(cellLeftAbove);
                }
            } else {
                if (rowSpaceScore > 0) {
                    currentCell.setScore(rowSpaceScore);
                    currentCell.setPrevCell(cellAbove);
                }
            }
        } else {
            if (matchOrMismatchScore >= colSpaceScore) {
                if (matchOrMismatchScore > 0) {
                    currentCell.setScore(matchOrMismatchScore);
                    currentCell.setPrevCell(cellLeftAbove);
                }
            } else {
                if (colSpaceScore > 0) {
                    currentCell.setScore(colSpaceScore);
                    currentCell.setPrevCell(cellToLeft);
                }
            }
        }
        if (currentCell.getScore() > highScoreCell.getScore()) {
            highScoreCell = currentCell;
        }
    }
    @Override
    protected Cell getInitialPointer(int row, int col) {
        return null;
    }
    @Override
    protected int getInitialScore(int row, int col) {
        return 0;
    }
}
4.编写测试类
public class test {
    public static void main(String[] args) {
        String s1 = "GCCCTAGCG";
        String s2 = "GCGCAATG";
        SmithWaterman sw = new SmithWaterman(s1, s2);
        sw.getScoreTable();
        String[] s = (String[]) sw.getTracebace();
        for (String str : s)
            System.out.println(str);
        NeedlemanWunsch nw = new NeedlemanWunsch(s1, s2);
        nw.getScoreTable();
        s = (String[]) sw.getTracebace();
        for (String str : s)
            System.out.println(str);
    }
}

1.两个算法通用的部分定义为一个抽象类，该类继承在LCS算法中定义的动态编程框架抽象类未进行优化，查找全局和局部比对的时间为 O(mn)，LCS及这两个算法都是动态编程的应用。动态编程还可用于矩阵链相乘、装配线规划和计算机象棋程序。使用动态编程能够解决的LCS及以上两个问题，这些问题都有共同的特征：每个问题的解都能用递归关系表示;用递归方法对这些递归关系的直接实现会造成解决方案效率低下，因为其中包含了对子问题的多次计算;一个问题的最优解能够用原始问题的子问题的最优解构造得到.

http://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm

http://en.wikipedia.org/wiki/Smith%E2%80%93Waterman_algorithm

http://www.w3school.com.cn/php/php_ref_mysql.asp

http://hi.baidu.com/%CF%FE%D4%C2%B7%C9%B7%C9/blog/item/693e1481c53417dd9123d93d.html