Pairwise Sequence Alignment (PSA) in java

Pairwise Sequence Alignment (PSA) 是一种在生物信息学中用于比较两个生物序列（例如DNA、RNA或蛋白质序列）的方法。它的目的是通过引入间隙（gap）来对齐序列，从而最大化序列之间的相似性。最常见的PSA方法包括全局对齐（Needleman-Wunsch算法）和局部对齐（Smith-Waterman算法）。

全局对齐（Global Alignment）

全局对齐旨在对齐整个序列，从第一个字符到最后一个字符。它适用于长度相近并且整体上相似的序列。

Needleman-Wunsch算法（全局对齐）

Needleman-Wunsch算法用于全局对齐，旨在对齐整个序列。它基于动态规划方法，并使用一个得分矩阵来计算两个序列的对齐得分。

算法步骤：

1. 初始化得分矩阵：

   • 创建一个二维得分矩阵 score，大小为 (m+1) x (n+1)，其中 m 和 n 分别是两个序列的长度。
   • 初始化矩阵的第一行和第一列，表示对齐空序列的代价。

2. 填充得分矩阵：

   • 根据匹配、错配和间隙得分填充矩阵。
   • 使用动态规划公式更新矩阵：score[i][j] = max(score[i-1][j-1] + match/mismatch, score[i-1][j] + gap, score[i][j-1] + gap)。

3. 回溯构建对齐结果：

   • 从矩阵的右下角开始回溯，构建两个序列的对齐结果。

Java实现Needleman-Wunsch算法

以下是Java实现全局对齐（Needleman-Wunsch算法）的代码：

public class NeedlemanWunsch {

    // 方法：对两个序列进行全局对齐
    public static void globalAlignment(String seq1, String seq2, int match, int mismatch, int gap) {
        int m = seq1.length();
        int n = seq2.length();

        // 创建并初始化得分矩阵
        int[][] score = new int[m + 1][n + 1];

        for (int i = 0; i <= m; i++) {
            score[i][0] = i * gap;
        }
        for (int j = 0; j <= n; j++) {
            score[0][j] = j * gap;
        }

        // 填充得分矩阵
        for (int i = 1; i <= m; i++) {
            for (int j = 1; j <= n; j++) {
                int matchMismatch = (seq1.charAt(i - 1) == seq2.charAt(j - 1)) ? match : mismatch;
                score[i][j] = Math.max(Math.max(score[i - 1][j - 1] + matchMismatch, score[i - 1][j] + gap), score[i][j - 1] + gap);
            }
        }

        // 输出对齐得分
        System.out.println("Alignment Score: " + score[m][n]);

        // 回溯构建对齐结果
        StringBuilder align1 = new StringBuilder();
        StringBuilder align2 = new StringBuilder();

        int i = m, j = n;
        while (i > 0 && j > 0) {
            if (score[i][j] == score[i - 1][j - 1] + ((seq1.charAt(i - 1) == seq2.charAt(j - 1)) ? match : mismatch)) {
                align1.append(seq1.charAt(i - 1));
                align2.append(seq2.charAt(j - 1));
                i--;
                j--;
            } else if (score[i][j] == score[i - 1][j] + gap) {
                align1.append(seq1.charAt(i - 1));
                align2.append('-');
                i--;
            } else {
                align1.append('-');
                align2.append(seq2.charAt(j - 1));
                j--;
            }
        }

        // 如果有剩余字符
        while (i > 0) {
            align1.append(seq1.charAt(i - 1));
            align2.append('-');
            i--;
        }
        while (j > 0) {
            align1.append('-');
            align2.append(seq2.charAt(j - 1));
            j--;
        }

        // 输出对齐结果
        System.out.println("Sequence 1: " + align1.reverse().toString());
        System.out.println("Sequence 2: " + align2.reverse().toString());
    }

    public static void main(String[] args) {
        String seq1 = "GATTACA";
        String seq2 = "GCATGCU";
        int match = 1;
        int mismatch = -1;
        int gap = -1;

        globalAlignment(seq1, seq2, match, mismatch, gap);
    }
}

代码说明

1. 初始化得分矩阵：

   • score[i][0] = i * gap：对齐第一个序列的前 i 个字符和空序列的得分。
   • score[0][j] = j * gap：对齐第二个序列的前 j 个字符和空序列的得分。

2. 填充得分矩阵：

   • 使用动态规划公式计算每个位置的得分。
   • matchMismatch 是根据字符是否相等决定的得分。

3. 回溯构建对齐结果：

   • 从右下角开始回溯，根据得分矩阵的值决定当前字符的对齐方式。
   • 当 i 或 j 大于 0 时，如果有剩余字符，需要将它们与空字符对齐。

4. 输出对齐结果：

   • 输出对齐得分和对齐后的两个序列。

示例输出

对于输入的序列 seq1 = "GATTACA" 和 seq2 = "GCATGCU"，代码将输出：

Alignment Score: 0
Sequence 1: GATTACA-
Sequence 2: GCAT-GCU

这段代码实现了全局序列对齐（Needleman-Wunsch算法），并能够输出对齐得分及对齐后的序列。

局部对齐（Local Alignment）

局部对齐旨在找到两个序列中最相似的片段，适用于长度差异较大或者只需要比较局部相似性的序列。

Smith-Waterman算法

这是实现局部对齐的经典算法，同样基于动态规划方法，但其填充和回溯规则与全局对齐不同。

关键特性：

1. 对齐局部片段：只对齐序列中最相似的局部区域。
2. 适用于长度差异大的序列：通常用于找到不同基因或蛋白质中的保守区域（conserved regions）。

Java实现Smith-Waterman算法：

public class SmithWaterman {

    // 方法：对两个序列进行局部对齐
    public static void localAlignment(String seq1, String seq2, int match, int mismatch, int gap) {
        int m = seq1.length();
        int n = seq2.length();

        // 创建并初始化得分矩阵
        int[][] score = new int[m + 1][n + 1];
        int maxScore = 0;
        int maxI = 0, maxJ = 0;

        // 填充得分矩阵
        for (int i = 1; i <= m; i++) {
            for (int j = 1; j <= n; j++) {
                int matchMismatch = (seq1.charAt(i - 1) == seq2.charAt(j - 1)) ? match : mismatch;
                score[i][j] = Math.max(0, Math.max(score[i - 1][j - 1] + matchMismatch, Math.max(score[i - 1][j] + gap, score[i][j - 1] + gap)));

                if (score[i][j] > maxScore) {
                    maxScore = score[i][j];
                    maxI = i;
                    maxJ = j;
                }
            }
        }

        // 输出局部对齐得分
        System.out.println("Local Alignment Score: " + maxScore);

        // 回溯构建对齐结果
        StringBuilder align1 = new StringBuilder();
        StringBuilder align2 = new StringBuilder();

        int i = maxI, j = maxJ;
        while (i > 0 && j > 0 && score[i][j] != 0) {
            if (score[i][j] == score[i - 1][j - 1] + ((seq1.charAt(i - 1) == seq2.charAt(j - 1)) ? match : mismatch)) {
                align1.append(seq1.charAt(i - 1));
                align2.append(seq2.charAt(j - 1));
                i--;
                j--;
            } else if (score[i][j] == score[i - 1][j] + gap) {
                align1.append(seq1.charAt(i - 1));
                align2.append('-');
                i--;
            } else {
                align1.append('-');
                align2.append(seq2.charAt(j - 1));
                j--;
            }
        }

        // 输出对齐结果
        System.out.println("Sequence 1: " + align1.reverse().toString());
        System.out.println("Sequence 2: " + align2.reverse().toString());
    }

    public static void main(String[] args) {
        String seq1 = "GATTACA";
        String seq2 = "GCATGCU";
        int match = 1;
        int mismatch = -1;
        int gap = -1;

        localAlignment(seq1, seq2, match, mismatch, gap);
    }
}

j

1. 初始化得分矩阵：

   • 不同于全局对齐，局部对齐的得分矩阵初始值为0，表示不对齐的得分。

2. 填充得分矩阵：

   • 使用动态规划公式填充矩阵，同时记录最大得分及其位置。

3. 回溯构建对齐结果：

   • 从最大得分位置开始回溯，直到遇到得分为0的位置。

4. 输出对齐结果：

   • 输出局部对齐得分和对齐后的两个序列。

示例输出

对于输入的序列 seq1 = "GATTACA" 和 seq2 = "GCATGCU"，代码将输出：

Local Alignment Score: 3
Sequence 1: ATT
Sequence 2: ATG

总结

• 全局对齐（Global Alignment）：用于对齐整个序列，适用于长度相近的序列，常用算法是Needleman-Wunsch算法。
• 局部对齐（Local Alignment）：用于找到最相似的局部区域，适用于长度差异较大的序列，常用算法是Smith-Waterman算法。