基因序列比较

目录

1.问题描述

2.解决问题所用的方法及基本思路

3.算法描述

4.算法时间空间复杂度分析

5.算法实例

6.代码

---

1.问题描述 

人类基因由4种核苷酸，分别用字母ACTG表示。要求编写一个程序，按以下规划比较两个基因序列并确定它们的相似程度。即给出两个基因序列AGTGATG和GTTAG，它们有多相似呢？测量两个基因的相似度一种方法称为对齐。使用对齐方法可以在基因的适当位置加入空格，让两个基因的长度相等，然后根据基因的分值矩阵计算分数。

比较AGTGATG与GTTAG

第一种对齐方案为：

首先可以给AGTATG插入一个空格得：AGTAT-G

         GTTAG插入3个空格即得：-GT--TAG

上面的匹配分值为:-3+5+5+(-2)+(-3)+5+(-3)+5=9.

第二种对齐方案为：

AGTGATG

-GTTA-G

得到的分值为：(-3)+5+5+(-2)+5+(-1)+5=14.

  当然还有其它对齐方式，但以上对齐方式是最优的，所以两个基因的相似度就为14。

2.解决问题所用的方法及基本思路 

 问题是判断基因序列的相似度，可以用动态规划和减治的方法进行求解。减治法每一次先求解出更小一规模的子问题，再将小规模问题的解进行合并。动态规划法是采用自下而上的方式求值，并把中间最优结果存储起来。

  根据题目的对齐方式，设两条基因序列分别为长度为len1的gene1序列和长度为len2的gene2序列，可以理解为要找一条由坐标（0,0）到坐标（len1-1，len2-1）的路径，规则是只能向下、向右或者向右下走，使得这条路径的值最高,方式如下图所示：

两序列对比的过程用二维数组scoreMatrix[i][j]显示，其中i表示的是第一个序列gene1的第i个元素，j表示的是第二个序列gene2的第j个元素，矩阵中的每一元素可表示第一个序列gene1已比对字母且第二个序列gene2已比对字母相同比对结果的最优者，右下格就是整个比对的最优结果。

  可以看到，原来的插值、匹配问题已经被转换成一个二维寻址问题，每一步的选择是向下、向右或者向右下走，使得这条路径的值最高。

  具体计算时每一格只用最优比对的得分表示，给出的长度为m的X序列和长度为n的Y序列，建立一个（m+2）*（n+2）的矩阵，函数scoreMatrix[i][j]递归求序列gene1[1...i]和gene2[1....j]的最佳比配的得分，建立以下状态转移方程：

scoreMatrix[i][j]=max{scoreMatrix[i-1][j-1]+Matchscore,scoreMatrix[i-1][j]+Gapscore(gene1,’-’),scoreMatrix[i][j-1]+Gapscore(‘-’,gene2)}

3.算法描述 

动态规划算法：calculateScoreMatrix(gene1, gene2, len1, len2)

//输入：基因序列gene1,gene2,基因序列的长度len1,len2，

//输出：分值矩阵，对齐路径及最大得分

初始化分值矩阵

For i←0 to m：

 For j←0 to n：

scoreMatrix[i][j]=0

For i←1 to m：

    scoreMatrix[i][0]=scoreMatrix[i-1][0]+Gapscoreindex(gene1[i-1])

For j←1 to n：

    scoreMatrix[0][j]=scoreMatrix[0][j-1]+Gapscoreindex(gene1[j-1])

For i←0 to m：

 For j←0 to n：

scoreMatrix[i][j]=max{scoreMatrix[i-1][j-1]+Matchscore,scoreMatrix[i-1][j]+Gapscore(gene1,’-’),scoreMatrix[i][j-1]+Gapscore(‘-’,gene2)}

减治算法：divideAndConquerAlignment(gene1,gene2,len1,len2)

// 输入：基因序列gene1,gene2,基因序列的长度len1,len2，

// 输出：左上部分的分值矩阵和对齐路径及右下部分的分值矩阵和对齐路径，最大得分

* scoreMatrix1 = calculateScoreMatrix(gene1, gene2, mid1, mid2)

For i←0 to mid1：

 For j←0 to mid2：

            printf(scoreMatrix1[i][j])

 traceback(gene1, gene2, mid1, mid2, scoreMatrix1)

* scoreMatrix2 = calculateScoreMatrix(gene1 + mid1, gene2 + mid2, len1 - mid1, len2 - mid2)

For i←0 to len1-mid1：

 For j←0 to len2-mid2：

       printf(scoreMatrix1[i][j])

 traceback(gene1 + mid1, gene2 + mid2, len1 - mid1, len2 - mid2, scoreMatrix2)

a←scoreMatrix1[mid1][mid2];

b←scoreMatrix2[len1-mid1][len2-mid2];

return (a+b)

4.算法时间空间复杂度分析 

基本运算：比较操作

基本运算的重复执行次数的表达式：C(n) = m*n

时间复杂度渐进表示：O(mn)

空间复杂度渐进表示：O(mn)

  当在序列中插入一个空字符串时，需要新建两个数组来存放新的序列，这步操作花费了m+n的时间。在初始化时有m+n次赋值操作，加上运算时m×n次赋值操作，一共(m×n)+m+n次赋值，所以最后的时间复杂度为O(mn)。

  因为运算时，创建了一个记录结果的表，所以空间复杂度为（m+2)×(n+2)，所以最后的空间复杂度为O(mn)。

5.算法实例 

6.代码 

 6.1 动态规划+回溯法

#include <stdio.h>
#include <stdlib.h>

int max(int a, int b, int c) {
    if (a >= b && a >= c) {
        return a;
    } else if (b >= a && b >= c) {
        return b;
    } else {
        return c;
    }
}

int max1(int a,int b)
{
    if(a>b)
        return a;
    else
        return b;
}

int min(int a,int b)
{
    if(a<b)
        return a;
    else
        return b;
}

int Gapscore(char* gene,int len)
{
    int i;
    int score = 0;
    for(i=0;i<len;i++)
    {
       if(gene[i]=='A')
           score += -3;
       else if(gene[i]=='C')
           score += -4;
       else if(gene[i]=='G')
           score += -2;
       else if(gene[i]=='T')
           score += -1;
    }
    return score;
}

int Gapscoreindex(char a)
{
    if(a=='A')
        return -3;
    else if(a=='C')
        return -4;
    else if(a=='G')
        return -2;
    else if(a=='T')
        return -1;
    else
        return 0;
}

int Matchscore(char a, char b)
{
    if (a == b) 
        return 5;
    else if ((a == 'A' && b == 'C') || (a == 'C' && b == 'A'))
        return -1;
    else if ((a == 'A' && b == 'G') || (a == 'G' && b == 'A'))
        return -2;
    else if ((a == 'A' && b == 'T') || (a == 'T' && b == 'A'))
        return -1;
    else if ((a == 'C' && b == 'G') || (a == 'G' && b == 'C'))
        return -3;
    else if ((a == 'C' && b == 'T') || (a == 'T' && b == 'C'))
        return -2;
    else if ((a == 'G' && b == 'T') || (a == 'T' && b == 'G'))
        return -2;
    else
        return 5;
}
void traceback(char* gene1, char* gene2, int len1, int len2, int scoreMatrix[len1+1][len2+1])
//初始化变量：定义两个字符数组align1和align2用于存储比对后的序列，初始化i和j为gene1和gene2的长度，初始化k为0。
    {
int i = len1;
    int j = len2;
    char align1[len1+len2+1];
    char align2[len1+len2+1];
    int k = 0;
// 进入循环：在循环中进行回溯操作，直到i和j都为0。
while (i > 0 || j > 0) {
//检查匹配情况：如果i和j都大于0，并且gene1[i-1]等于gene2[j-1]，则说明当前位置的字符匹配，将其添加到align1和align2中，并将i和j分别减1，k加1。
        if (i > 0 && j > 0 && gene1[i-1] == gene2[j-1])
 {
            align1[k] = gene1[i-1];
            align2[k] = gene2[j-1];
            i--;
            j--;
            k++;
        } //检查得分情况：如果i和j都大于0，并且当前位置的得分等于scoreMatrix[i-1][j-1] + matchScore，说明当前位置的得分是通过匹配得到的，将对应的字符添加到align1和align2中，并将i和j分别减1，k加1。
        else if (i > 0 && j > 0)
 {
            int matchScore = Matchscore(gene1[i-1], gene2[j-1]);
            int gapScore1 = Gapscore(&gene1[i-1], 1);
            int gapScore2 = Gapscore(&gene2[j-1], 1);
            if (scoreMatrix[i][j] == scoreMatrix[i-1][j-1] + matchScore) {
                align1[k] = gene1[i-1];
                align2[k] = gene2[j-1];
                i--;
                j--;
                k++;
            }
//检查gap情况：如果i大于0，并且当前位置的得分等于scoreMatrix[i-1][j] + gapScore1，说明当前位置的得分是通过在gene1中插入一个gap得到的，将gene1中的字符添加到align1中，将'-'添加到align2中，并将i减1，k加1。
 else if (scoreMatrix[i][j] == scoreMatrix[i-1][j] + gapScore1) 
{
                align1[k] = gene1[i-1];
                align2[k] = '-';
                i--;
                k++;
            } 
//检查gap情况：如果j大于0，并且当前位置的得分等于scoreMatrix[i][j-1] + gapScore2，说明当前位置的得分是通过在gene2中插入一个gap得到的，将gene2中的字符添加到align2中，将'-'添加到align1中，并将j减1，k加1。
else if (scoreMatrix[i][j] == scoreMatrix[i][j-1] + gapScore2) {
                align1[k] = '-';
                align2[k] = gene2[j-1];
                j--;
                k++;
            }
        } 
// 处理剩余字符：如果i大于0，说明gene1还有剩余字符，将剩余字符添加到align1中，将'-'添加到align2中，并将i减1，k加1。
else if (i > 0)
 {
            int gapScore1 = Gapscore(&gene1[i-1], 1);
            align1[k] = gene1[i-1];
            align2[k] = '-';
            i--;
            k++;
        } 
//处理剩余字符：如果j大于0，说明gene2还有剩余字符，将剩余字符添加到align2中，将'-'添加到align1中，并将j减1，k加1。
else
 {
            int gapScore2 = Gapscore(&gene2[j-1], 1);
            align1[k] = '-';
            align2[k] = gene2[j-1];
            j--;
            k++;
        }
    }

    // 反转字符串,将align1和align2中的字符反转，使其按照正常顺序排列。

    int m;
for (m = 0; m < k/2; m++)
 {
        char temp = align1[m];
        align1[m] = align1[k-m-1];
        align1[k-m-1] = temp;

        temp = align2[m];
        align2[m] = align2[k-m-1];
        align2[k-m-1] = temp;
    }

    // 输出对齐结果
    printf("%s\n%s\n", align1, align2);
}

int** calculateScoreMatrix(char* gene1, char* gene2, int len1, int len2) {
    // 创建分值矩阵
    int** scoreMatrix = (int**)malloc((len1 + 1) * sizeof(int*));
    int i, j;

    for (i = 0; i <= len1; i++) {
        scoreMatrix[i] = (int*)malloc((len2 + 1) * sizeof(int));
    }

    // 初始化分值矩阵
    for (i = 0; i <= len1; i++) {
        for (j = 0; j <= len2; j++) {
            scoreMatrix[i][j] = 0;
        }
    }

    // 计算第一列的分值
    for (i = 1; i <= len1; i++) {
        int gapScore1 = Gapscoreindex(gene1[i - 1]);
        scoreMatrix[i][0] = scoreMatrix[i - 1][0] + gapScore1;
    }

    // 计算第一行的分值
    for (j = 1; j <= len2; j++) {
        int gapScore2 = Gapscoreindex(gene2[j - 1]);
        scoreMatrix[0][j] = scoreMatrix[0][j - 1] + gapScore2;
    }

    // 计算其它位置的分值
    for (i = 1; i <= len1; i++) {
        for (j = 1; j <= len2; j++) {
            int matchScore = Matchscore(gene1[i - 1], gene2[j - 1]);
            int gapScore1 = Gapscore(&gene1[i - 1], 1);
            int gapScore2 = Gapscore(&gene2[j - 1], 1);
            scoreMatrix[i][j] = max(scoreMatrix[i - 1][j - 1] + matchScore,
                                    scoreMatrix[i - 1][j] + gapScore1,
                                    scoreMatrix[i][j - 1] + gapScore2);
        }
    }
    printf("Score Matrix:\n");
for (i = 0; i <= len1; i++)
 {
        for (j = 0; j <= len2; j++)
 {
            printf("%d ", scoreMatrix[i][j]);
         }
        printf("\n");
    }
   
}
int main() {
    char gene1[] = "AGTGATG";
    char gene2[] = "GTTAG";
    int len1 = sizeof(gene1) - 1;
    int len2 = sizeof(gene2) - 1;
    int matchScore, gapScore1, gapScore2;
    int scoreMatrix[len1 + 1][len2 + 1];
    int i, j;

 calculateScoreMatrix(gene1, gene2, len1, len2);
    
   

    return 0;
}

6.2 减治法 

#include<stdio.h>
#include<stdlib.h>
int max(int a, int b, int c) 
{
if (a >= b && a >= c) 
{
        return a;
} 
else if (b >= a && b >= c)
 {
        return b;
} 
else 
{
        return c;
    }
}

int max1(int a,int b)
{
    if(a>b)
        return a;
    else
        return b;
}

int min(int a,int b)
{
    if(a<b)
        return a;
    else
        return b;
}

int Gapscore(char* gene,int len)
{
    int i;
    int score = 0;
    for(i=0;i<len;i++)
    {
       if(gene[i]=='A')
           score += -3;
       else if(gene[i]=='C')
           score += -4;
       else if(gene[i]=='G')
           score += -2;
       else if(gene[i]=='T')
           score += -1;
    }
    return score;
}

int Gapscoreindex(char a)
{
    if(a=='A')
        return -3;
    else if(a=='C')
        return -4;
    else if(a=='G')
        return -2;
    else if(a=='T')
        return -1;
    else
        return 0;
}

int Matchscore(char a, char b)
{
    if (a == b) 
        return 5;
    else if ((a == 'A' && b == 'C') || (a == 'C' && b == 'A'))
        return -1;
    else if ((a == 'A' && b == 'G') || (a == 'G' && b == 'A'))
        return -2;
    else if ((a == 'A' && b == 'T') || (a == 'T' && b == 'A'))
        return -1;
    else if ((a == 'C' && b == 'G') || (a == 'G' && b == 'C'))
        return -3;
    else if ((a == 'C' && b == 'T') || (a == 'T' && b == 'C'))
        return -2;
    else if ((a == 'G' && b == 'T') || (a == 'T' && b == 'G'))
        return -2;
    else
        return 5;
}
int** calculateScoreMatrix(char* gene1, char* gene2, int len1, int len2) {
    // 创建分值矩阵
    int** scoreMatrix = (int**)malloc((len1 + 1) * sizeof(int*));
    int i, j;

for (i = 0; i <= len1; i++)
 {
        scoreMatrix[i] = (int*)malloc((len2 + 1) * sizeof(int));
     }

    // 初始化分值矩阵
for (i = 0; i <= len1; i++)
 {
        for (j = 0; j <= len2; j++)
 {
            scoreMatrix[i][j] = 0;
         }
    }

    // 计算第一列的分值
for (i = 1; i <= len1; i++)
 {
        int gapScore1 = Gapscoreindex(gene1[i - 1]);
        scoreMatrix[i][0] = scoreMatrix[i - 1][0] + gapScore1;
     }

    // 计算第一行的分值
for (j = 1; j <= len2; j++)
 {
        int gapScore2 = Gapscoreindex(gene2[j - 1]);
        scoreMatrix[0][j] = scoreMatrix[0][j - 1] + gapScore2;
     }

    // 计算其它位置的分值
for (i = 1; i <= len1; i++)
 {
        for (j = 1; j <= len2; j++)
 {
            int matchScore = Matchscore(gene1[i - 1], gene2[j - 1]);
            int gapScore1 = Gapscore(&gene1[i - 1], 1);
            int gapScore2 = Gapscore(&gene2[j - 1], 1);
            scoreMatrix[i][j] = max(scoreMatrix[i - 1][j - 1] + matchScore, scoreMatrix[i - 1][j] + gapScore1, scoreMatrix[i][j - 1] + gapScore2);
        }
    }
    return scoreMatrix;
}
void traceback(char* gene1, char* gene2, int len1, int len2, int** scoreMatrix) 
{
//初始化变量：定义两个字符数组align1和align2用于存储比对后的序列，初始化i和j为gene1和gene2的长度，初始化k为0。
    {
int i = len1;
    int j = len2;
    char align1[len1+len2+1];
    char align2[len1+len2+1];
    int k = 0;
// 进入循环：在循环中进行回溯操作，直到i和j都为0。
while (i > 0 || j > 0) {
//检查匹配情况：如果i和j都大于0，并且gene1[i-1]等于gene2[j-1]，则说明当前位置的字符匹配，将其添加到align1和align2中，并将i和j分别减1，k加1。
        if (i > 0 && j > 0 && gene1[i-1] == gene2[j-1])
 {
            align1[k] = gene1[i-1];
            align2[k] = gene2[j-1];
            i--;
            j--;
            k++;
        } //检查得分情况：如果i和j都大于0，并且当前位置的得分等于scoreMatrix[i-1][j-1] + matchScore，说明当前位置的得分是通过匹配得到的，将对应的字符添加到align1和align2中，并将i和j分别减1，k加1。
        else if (i > 0 && j > 0)
 {
            int matchScore = Matchscore(gene1[i-1], gene2[j-1]);
            int gapScore1 = Gapscore(&gene1[i-1], 1);
            int gapScore2 = Gapscore(&gene2[j-1], 1);
            if (scoreMatrix[i][j] == scoreMatrix[i-1][j-1] + matchScore) {
                align1[k] = gene1[i-1];
                align2[k] = gene2[j-1];
                i--;
                j--;
                k++;
            }
//检查gap情况：如果i大于0，并且当前位置的得分等于scoreMatrix[i-1][j] + gapScore1，说明当前位置的得分是通过在gene1中插入一个gap得到的，将gene1中的字符添加到align1中，将'-'添加到align2中，并将i减1，k加1。
 else if (scoreMatrix[i][j] == scoreMatrix[i-1][j] + gapScore1) 
{
                align1[k] = gene1[i-1];
                align2[k] = '-';
                i--;
                k++;
            } 
//检查gap情况：如果j大于0，并且当前位置的得分等于scoreMatrix[i][j-1] + gapScore2，说明当前位置的得分是通过在gene2中插入一个gap得到的，将gene2中的字符添加到align2中，将'-'添加到align1中，并将j减1，k加1。
else if (scoreMatrix[i][j] == scoreMatrix[i][j-1] + gapScore2) {
                align1[k] = '-';
                align2[k] = gene2[j-1];
                j--;
                k++;
            }
        } 
// 处理剩余字符：如果i大于0，说明gene1还有剩余字符，将剩余字符添加到align1中，将'-'添加到align2中，并将i减1，k加1。
else if (i > 0)
 {
            int gapScore1 = Gapscore(&gene1[i-1], 1);
            align1[k] = gene1[i-1];
            align2[k] = '-';
            i--;
            k++;
        } 
//处理剩余字符：如果j大于0，说明gene2还有剩余字符，将剩余字符添加到align2中，将'-'添加到align1中，并将j减1，k加1。
else
 {
            int gapScore2 = Gapscore(&gene2[j-1], 1);
            align1[k] = '-';
            align2[k] = gene2[j-1];
            j--;
            k++;
        }
    }

    // 反转字符串,将align1和align2中的字符反转，使其按照正常顺序排列。

    int m;
for (m = 0; m < k/2; m++)
 {
        char temp = align1[m];
        align1[m] = align1[k-m-1];
        align1[k-m-1] = temp;

        temp = align2[m];
        align2[m] = align2[k-m-1];
        align2[k-m-1] = temp;
    }

    // 输出对齐结果
    printf("%s\n%s\n", align1, align2);
}

   

int divideAndConquerAlignment(char* gene1, char* gene2, int len1, int len2) 
{
if (len1 == 0 || len2 == 0) 
{
        // 基本情况：一个或两个基因序列为空
        printf("Alignment:\n%s\n%s\n", gene1, gene2);
        printf("Max Score: 0\n");
        return 0;
    }

if (len1 == 1 && len2 == 1)
 {
        // 基本情况：两个基因序列只有一个字符
        int matchScore = Matchscore(gene1[0], gene2[0]);
        printf("Alignment:\n%c\n%c\n", gene1[0], gene2[0]);
        printf("Max Score: %d\n", matchScore);
        return 0;
    }

    // 分解问题为更小的子问题
    int mid1 = len1 / 2;
    int mid2 = len2 / 2;

    // 计算左上部分的分值矩阵和对齐路径
    int** scoreMatrix1 = calculateScoreMatrix(gene1, gene2, mid1, mid2);
    printf("Score Matrix (Left Top):\n");
for (int i = 0; i <= mid1; i++)
 {
        for (int j = 0; j <= mid2; j++)
 {
            printf("%d ", scoreMatrix1[i][j]);
         }
         printf("\n");
    }
    printf("Alignment Path (Left Top):\n");
    traceback(gene1, gene2, mid1, mid2, scoreMatrix1);
    printf("\n");

    // 计算右下部分的分值矩阵和对齐路径
    int** scoreMatrix2 = calculateScoreMatrix(gene1 + mid1, gene2 + mid2, len1 - mid1, len2 - mid2);
    printf("Score Matrix (Right Bottom):\n");
for (int i = 0; i <= len1 - mid1; i++)
 {
        for (int j = 0; j <= len2 - mid2; j++)
 {
            printf("%d ", scoreMatrix2[i][j]);
        }
        printf("\n");
    }
    printf("Alignment Path (Right Bottom):\n");
    traceback(gene1 + mid1, gene2 + mid2, len1 - mid1, len2 - mid2, scoreMatrix2);
    printf("\n");
    int a=scoreMatrix1[mid1][mid2];
    int b=scoreMatrix2[len1-mid1][len2-mid2];
    // 释放分值矩阵的内存
for (int i = 0; i <= mid1; i++) 
{
        free(scoreMatrix1[i]);
    }
    free(scoreMatrix1);

for (int i = 0; i <= len1 - mid1; i++)
 {
        free(scoreMatrix2[i]);
     }
    free(scoreMatrix2);
    return (a+b);
}

int main()
 {
    char gene1[] = "AGTGATG";
    char gene2[] = "GTTAG";
    int len1 = sizeof(gene1) - 1;
    int len2 = sizeof(gene2) - 1;
    int m=divideAndConquerAlignment(gene1, gene2, len1, len2);
    printf("Max Score: %d\n", m);
    return 0;
}