天天看點

DNA序列對齊問題

一、問題描述

該問題在算法導論中引申自求解兩個DNA序列相似度的問題。

可以從很多角度定義兩個DNA序列的相似度,其中有一種定義方法就是通過序列對齊的方式來定義其相似度。

給定兩個DNA序列A和B,對齊的方式是将空格分别插入到A和B序列中,得到具有相同長度的對齊後的序列C和D;空格可以插入到任意的位置(包括兩端),但是相同位置不能同時為空格,也即是不存在C[i]和D[i]同時為空格的情況。然後為對齊後的序列的每個位置打分,總分為每個位置得分之和,具體的打分規則如下:

a、如果C[i] == D[i]且都不是空格,得3分;

b、如果C[i] != D[j]且都不是空格,得1分;

c、如果C[i] 或者D[i]是空格,得0分。

求給定原序列A和B的一個對齊方案,使得該對齊方案的總分最高。

例如,序列原序列A和B如下:

String strA = "GATC";
        String strB = "ATCG";      

則其中一個對齊方案如下:

GATC*
*ATCG      

該方案總得分score=2*0+3*3 = 9分。

實際上這是最優的對齊方案,在所有的對齊方案中總得分最高為9分。

二、問題分析

為了用更加簡單的方式來表示對齊的方案,我們嘗試用一些特定的字元記号來表示對齊方案,對此,首先做一個約定,對于打分規則:

1、情況a用“=”字元标記;

2、情況b用“~”字元标記;

3、情況c用“*”字元标記,但是情況c實際上可以細分為兩種情況:C[i]為空格時用“+”标記,D[i]為空格時用“-”号标記。這樣用“+”和“-”細分的表示相比于統一用“*”來表示,本質的差別在于讓對齊方案具有所謂的“方向性”,後面會看到這樣的細分對于算法的實作有一定的好處。

有了這樣的約定,可以将一個對齊方案用這些字元表示出來,該字元串稱之為一個對齊規則字元串R。

例如上面的例子中,對齊規則就可以用字元串“-===+”來表示。

可以推斷,任何兩個原序列的對齊規則字元串R的長度必然滿足:

DNA序列對齊問題

隻要能夠求得最優對齊方案的對齊規則字元串,就可以計算出最高分數,還可以還原出各自的對齊序列。

考察該問題的最優子結構性質,與最長公共子序列思考的角度比較類似,

用C(i,j)表示序列A[0]...A[i]和序列B[0]...B[j]的最優對齊方案的得分,不難得出其初始條件和遞推求解式:

DNA序列對齊問題

用R(i,j)表示序列A[0]...A[i]和序列B[0]...B[j]的最優對齊方案的對齊規則字元串,結合上面的遞推求解式,不難推出對齊規則字元串的運算規則:

DNA序列對齊問題

三、算法實作

package agdp;
public class Alignment {
    //根據對齊骨規則生成相應的對齊後的字元串
    private static String[] generate(String base,String ...origin){
        int num = origin.length;
        String[] align = new String[num];
        for (int i = 0; i < num; i++) {
            if (origin[i].length() == base.length()) {
                align[i] = origin[i];
            }else {//base.length()隻能是等于或者大于兩個原字元串的長度
                String tmp = "";
                for (int j = 0,k = 0; j < base.length(); j++) {
                    if (base.charAt(j) == '+') {
                        if (i == 0) {tmp = tmp+"*";}
                        else{tmp = tmp+origin[i].charAt(k++);}
                    }else if(base.charAt(j) == '-'){
                        if (i == 0) {tmp = tmp+origin[i].charAt(k++);}
                        else{tmp = tmp+"*";}
                    }
                    else {
                        tmp = tmp+origin[i].charAt(k++);
                    }
                }
                align[i] = tmp;
            }
        }
        return align;
    }
    public static String align(String strA,String strB){
        int m = strA.length(),n = strB.length(),tmp;
        //aux數組記錄子問題的最有對齊方案的分數,也即子問題的最高分數。
        int[][] aux = new int[m+1][n+1];
        //rule數組記錄對齊方案,分别用"+"、"-"、"="和"~"記錄四種情況。
        String[][] rule = new String[m+1][n+1];
        //rule初始化
        rule [0][0] = "";
        for (int i = 1; i < m+1; i++) {
            rule[i][0] = rule[i-1][0]+"-";
        }
        for (int i = 1; i < n+1; i++) {
            rule[0][i] = rule[0][i-1]+"+";
        }
        for (int i = 1; i <= m; i++) {
            for (int j = 1; j <= n; j++) {
                if (strA.charAt(i-1) == strB.charAt(j-1)) {
                    aux[i][j] = aux[i-1][j-1]+3;
                    rule[i][j] = rule[i-1][j-1] + "=";//A[i]==B[j]:->"="
                }else {
                    tmp = Math.max(Math.max(aux[i-1][j], aux[i][j-1]), aux[i-1][j-1]+1);
                    aux[i][j] = tmp;
                    if (tmp == aux[i-1][j-1]-1) {//A[i]!=B[j]且A[i]和    B[j]不為空字元:->"~"
                        rule[i][j] = rule[i-1][j-1]+"~";
                    }else if(tmp == aux[i-1][j]-2){//B[i]為空字元:->"-"
                        rule[i][j] = rule[i-1][j]+"-";
                    }else{
                        rule[i][j] = rule[i][j-1]+"+";//A[i]為空字元:->"+"
                    }
                }
            }
        }
        //格式化輸出aux數組
        for (int i = 0; i < m+1; i++) {
            for (int j = 0; j < n+1; j++) {
                System.out.format("%3d",aux[i][j]);
            }
            System.out.println();
        }
        //格式化輸出rule數組
        for (int i = 0; i < m+1; i++) {
            for (int j = 0; j < n+1; j++) {
                System.out.format("%-15s",rule[i][j]);
            }
            System.out.println();
        }
        //傳回最優的對齊方法對應的規則
        return rule[m][n];
    }
    //根據規則字元串計算分數
    public static int getScore(String ruleStr){
        int score = 0;
        for (int i = 0; i < ruleStr.length(); i++) {
            if (ruleStr.charAt(i) == '=') {
                score += 3;
            }else if (ruleStr.charAt(i) == '~') {
                score += 1;
            }
        }
        return score;
    }
    public static void main(String[] args) {
        // TODO Auto-generated method stub
        int[] scoreAry = {3,1,0,0};
//        String strA = "GATCGGCAT";
//        String strB = "CAATGTGAATC";
        String strA = "GATC";
        String strB = "ATCG";
//        String strA = "GAC";
//        String strB = "ATCG";
        String ruleStr = align(strA, strB);
        System.out.println(ruleStr);
        int score = getScore(ruleStr);
        System.out.println(score);
        String[] alignStr = generate(ruleStr, strA,strB); 
        for(String str:alignStr){
            System.out.println(str);
        }
    }
}      

還是以原始序列“GATC”和“ATCG”為例:

其子問題的得分的計算如下:

DNA序列對齊問題

子問題的對齊規則字元串的計算如下:

DNA序列對齊問題

需要特别注意的是,用“+”和“-”号來區分打分情況c後,對齊規則字元串是具有“方向性”的,也就是說對齊規則“-===+”是指從A->B方向的對齊規則。那如果需要B->A的對齊規則,隻需要将對齊規則的字元串中+“和”-”互相替換即可。

實際上,從DNA序列對齊問題過渡到編輯距離問題是很比較自然的。本文也有意識的将這兩個問題聯系在一起,編輯距離問題見下一篇博文。

參考資料:

算法導論.第十五章 習題15-5

轉載請注明原文出處:

http://www.cnblogs.com/qcblog/p/7820140.html

繼續閱讀