Problem Statement

Genes consist of four nucleotide bases denoted by A, C, T, and G. The task involves comparing two gene sequences and determining their similarity by finding an optimal alignment. The alignment process inserts gaps into the sequences to make them equal in length, then computes a score based on a scoring matrix.

For sequences AGTGATG and GTTAG, we can illustrate two possible alignments:

Alignment Option 1:

AGTAT-G
-GT--TAG

Score: -3+5+5+(-2)+(-3)+5+(-3)+5 = 9

Alignment Option 2:

AGTGATG
-GTTA-G

Score: (-3)+5+5+(-2)+5+(-1)+5 = 14

The optimal alignment yields a similarity score of 14.

Approach and Methodology

The problem of measuring sequence similarity reduces to finding a maximum-scoring path through a 2D grid. Let gene1 and gene2 be sequences of lengths len1 and len2 respectively. Each cell in the matrix represents the best alignment score between prefixes of the sequences.

A two-dimensional array scoreMatrix[i][j] stores the optimal alignment score for gene1[1...i] versus gene2[1...j]. Each step in the DP table involves choosing between: diagonal movement (match/mismatch), downward movement (gap in sequence 2), or rightward movement (gap in sequence 1).

The recurrence relation becomes:

scoreMatrix[i][j] = max(
    scoreMatrix[i-1][j-1] + matchScore(gene1[i], gene2[j]),
    scoreMatrix[i-1][j] + gapPenalty(gene1[i], '-'),
    scoreMatrix[i][j-1] + gapPenalty('-', gene2[j])
)

Algorithm Specification

Dynamic Programming Approach

calculateScoreMatrix(seq1, seq2, m, n)
    create matrix dp[m+1][n+1]
    
    for i from 0 to m:
        for j from 0 to n:
            dp[i][j] = 0
    
    // initialize first column with gap penalties
    for i from 1 to m:
        dp[i][0] = dp[i-1][0] + gapCost(seq1[i-1])
    
    // initialize first row with gap penalties  
    for j from 1 to n:
        dp[0][j] = dp[0][j-1] + gapCost(seq2[j-1])
    
    // fill remaining cells
    for i from 1 to m:
        for j from 1 to n:
            match = dp[i-1][j-1] + matchScore(seq1[i-1], seq2[j-1])
            del = dp[i-1][j] + gapCost(seq1[i-1])
            ins = dp[i][j-1] + gapCost(seq2[j-1])
            dp[i][j] = max(match, del, ins)
    
    return dp

Divide-and-Conquer Approach

divideAndConquer(seq1, seq2, len1, len2)
    if len1 == 0 or len2 == 0:
        return 0
    
    if len1 == 1 and len2 == 1:
        return matchScore(seq1[0], seq2[0])
    
    mid1 = len1 / 2
    mid2 = len2 / 2
    
    // process upper-left quadrant
    dp1 = calculateScoreMatrix(seq1, seq2, mid1, mid2)
    printMatrix(dp1, mid1+1, mid2+1)
    traceback(seq1, seq2, mid1, mid2, dp1)
    
    // process lower-right quadrant  
    dp2 = calculateScoreMatrix(seq1+mid1, seq2+mid2, len1-mid1, len2-mid2)
    printMatrix(dp2, len1-mid1+1, len2-mid2+1)
    traceback(seq1+mid1, seq2+mid2, len1-mid1, len2-mid2, dp2)
    
    partA = dp1[mid1][mid2]
    partB = dp2[len1-mid1][len2-mid2]
    
    free(dp1); free(dp2)
    return partA + partB

Complexity Analysis

The algorithm performs O(m×n) comparisons, where m and n represent the lengths of the input sequences. Each cell in the DP table requires constant time to evaluate, resulting in quadratic time complexity O(mn).

Space requirements stem from storing the complete DP table, which occupies (m+1)×(n+1) integer entries. This yields O(mn) space complexity.

Initialization of boundary rows contributes O(m+n) operations, but this remains dominated by the main computational work of O(mn).

Implementation

Dynamic Programming with Backtracking

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

int maximum(int a, int b, int c) {
    int result = a;
    if (b > result) result = b;
    if (c > result) result = c;
    return result;
}

int gapPenalty(char nucleotide) {
    switch(nucleotide) {
        case 'A': return -3;
        case 'C': return -4;
        case 'G': return -2;
        case 'T': return -1;
        default: return 0;
    }
}

int substitutionScore(char a, char b) {
    if (a == b) return 5;
    
    if ((a == 'A' && b == 'C') || (a == 'C' && b == 'A')) return -1;
    if ((a == 'A' && b == 'G') || (a == 'G' && b == 'A')) return -2;
    if ((a == 'A' && b == 'T') || (a == 'T' && b == 'A')) return -1;
    if ((a == 'C' && b == 'G') || (a == 'G' && b == 'C')) return -3;
    if ((a == 'C' && b == 'T') || (a == 'T' && b == 'C')) return -2;
    if ((a == 'G' && b == 'T') || (a == 'T' && b == 'G')) return -2;
    
    return 5;
}

int** buildDPTable(char* s1, char* s2, int rows, int cols) {
    int** table = malloc((rows + 1) * sizeof(int*));
    for (int i = 0; i <= rows; i++)
        table[i] = malloc((cols + 1) * sizeof(int));
    
    for (int i = 0; i <= rows; i++)
        for (int j = 0; j <= cols; j++)
            table[i][j] = 0;
    
    for (int i = 1; i <= rows; i++)
        table[i][0] = table[i-1][0] + gapPenalty(s1[i-1]);
    
    for (int j = 1; j <= cols; j++)
        table[0][j] = table[0][j-1] + gapPenalty(s2[j-1]);
    
    for (int i = 1; i <= rows; i++) {
        for (int j = 1; j <= cols; j++) {
            int diag = table[i-1][j-1] + substitutionScore(s1[i-1], s2[j-1]);
            int up = table[i-1][j] + gapPenalty(s1[i-1]);
            int left = table[i][j-1] + gapPenalty(s2[j-1]);
            table[i][j] = maximum(diag, up, left);
        }
    }
    
    return table;
}

void backtrackAlignment(char* s1, char* s2, int rows, int cols, int** dp) {
    char result1[rows + cols + 1];
    char result2[rows + cols + 1];
    int idx = 0;
    int i = rows, j = cols;
    
    while (i > 0 || j > 0) {
        if (i > 0 && j > 0 && s1[i-1] == s2[j-1]) {
            result1[idx] = s1[i-1];
            result2[idx] = s2[j-1];
            i--; j--; idx++;
        } else if (i > 0 && j > 0) {
            int match = substitutionScore(s1[i-1], s2[j-1]);
            int gapUp = gapPenalty(s1[i-1]);
            int gapLeft = gapPenalty(s2[j-1]);
            
            if (dp[i][j] == dp[i-1][j-1] + match) {
                result1[idx] = s1[i-1];
                result2[idx] = s2[j-1];
                i--; j--; idx++;
            } else if (dp[i][j] == dp[i-1][j] + gapUp) {
                result1[idx] = s1[i-1];
                result2[idx] = '-';
                i--; idx++;
            } else {
                result1[idx] = '-';
                result2[idx] = s2[j-1];
                j--; idx++;
            }
        } else if (i > 0) {
            result1[idx] = s1[i-1];
            result2[idx] = '-';
            i--; idx++;
        } else {
            result1[idx] = '-';
            result2[idx] = s2[j-1];
            j--; idx++;
        }
    }
    
    for (int k = 0; k < idx / 2; k++) {
        char tmp = result1[k];
        result1[k] = result1[idx - k - 1];
        result1[idx - k - 1] = tmp;
        
        tmp = result2[k];
        result2[k] = result2[idx - k - 1];
        result2[idx - k - 1] = tmp;
    }
    
    result1[idx] = '\0';
    result2[idx] = '\0';
    printf("%s\n%s\n", result1, result2);
}

int main() {
    char seq1[] = "AGTGATG";
    char seq2[] = "GTTAG";
    int len1 = strlen(seq1);
    int len2 = strlen(seq2);
    
    int** dpTable = buildDPTable(seq1, seq2, len1, len2);
    
    printf("DP Matrix:\n");
    for (int i = 0; i <= len1; i++) {
        for (int j = 0; j <= len2; j++)
            printf("%d ", dpTable[i][j]);
        printf("\n");
    }
    
    printf("\nOptimal Alignment:\n");
    backtrackAlignment(seq1, seq2, len1, len2, dpTable);
    printf("\nFinal Score: %d\n", dpTable[len1][len2]);
    
    for (int i = 0; i <= len1; i++)
        free(dpTable[i]);
    free(dpTable);
    
    return 0;
}

Divide-and-Conquer Implementation

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

int maximum(int a, int b, int c) {
    int result = a;
    if (b > result) result = b;
    if (c > result) result = c;
    return result;
}

int gapPenalty(char nucleotide) {
    switch(nucleotide) {
        case 'A': return -3;
        case 'C': return -4;
        case 'G': return -2;
        case 'T': return -1;
        default: return 0;
    }
}

int substitutionScore(char a, char b) {
    if (a == b) return 5;
    
    if ((a == 'A' && b == 'C') || (a == 'C' && b == 'A')) return -1;
    if ((a == 'A' && b == 'G') || (a == 'G' && b == 'A')) return -2;
    if ((a == 'A' && b == 'T') || (a == 'T' && b == 'A')) return -1;
    if ((a == 'C' && b == 'G') || (a == 'G' && b == 'C')) return -3;
    if ((a == 'C' && b == 'T') || (a == 'T' && b == 'C')) return -2;
    if ((a == 'G' && b == 'T') || (a == 'T' && b == 'G')) return -2;
    
    return 5;
}

int** buildDPTable(char* s1, char* s2, int rows, int cols) {
    int** table = malloc((rows + 1) * sizeof(int*));
    for (int i = 0; i <= rows; i++)
        table[i] = malloc((cols + 1) * sizeof(int));
    
    for (int i = 0; i <= rows; i++)
        for (int j = 0; j <= cols; j++)
            table[i][j] = 0;
    
    for (int i = 1; i <= rows; i++)
        table[i][0] = table[i-1][0] + gapPenalty(s1[i-1]);
    
    for (int j = 1; j <= cols; j++)
        table[0][j] = table[0][j-1] + gapPenalty(s2[j-1]);
    
    for (int i = 1; i <= rows; i++) {
        for (int j = 1; j <= cols; j++) {
            int diag = table[i-1][j-1] + substitutionScore(s1[i-1], s2[j-1]);
            int up = table[i-1][j] + gapPenalty(s1[i-1]);
            int left = table[i][j-1] + gapPenalty(s2[j-1]);
            table[i][j] = maximum(diag, up, left);
        }
    }
    
    return table;
}

void backtrackAlignment(char* s1, char* s2, int rows, int cols, int** dp) {
    char result1[rows + cols + 1];
    char result2[rows + cols + 1];
    int idx = 0;
    int i = rows, j = cols;
    
    while (i > 0 || j > 0) {
        if (i > 0 && j > 0 && s1[i-1] == s2[j-1]) {
            result1[idx] = s1[i-1];
            result2[idx] = s2[j-1];
            i--; j--; idx++;
        } else if (i > 0 && j > 0) {
            int match = substitutionScore(s1[i-1], s2[j-1]);
            int gapUp = gapPenalty(s1[i-1]);
            int gapLeft = gapPenalty(s2[j-1]);
            
            if (dp[i][j] == dp[i-1][j-1] + match) {
                result1[idx] = s1[i-1];
                result2[idx] = s2[j-1];
                i--; j--; idx++;
            } else if (dp[i][j] == dp[i-1][j] + gapUp) {
                result1[idx] = s1[i-1];
                result2[idx] = '-';
                i--; idx++;
            } else {
                result1[idx] = '-';
                result2[idx] = s2[j-1];
                j--; idx++;
            }
        } else if (i > 0) {
            result1[idx] = s1[i-1];
            result2[idx] = '-';
            i--; idx++;
        } else {
            result1[idx] = '-';
            result2[idx] = s2[j-1];
            j--; idx++;
        }
    }
    
    for (int k = 0; k < idx / 2; k++) {
        char tmp = result1[k];
        result1[k] = result1[idx - k - 1];
        result1[idx - k - 1] = tmp;
        
        tmp = result2[k];
        result2[k] = result2[idx - k - 1];
        result2[idx - k - 1] = tmp;
    }
    
    result1[idx] = '\0';
    result2[idx] = '\0';
    printf("%s\n%s\n", result1, result2);
}

void printMatrix(int** mat, int rows, int cols) {
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++)
            printf("%d ", mat[i][j]);
        printf("\n");
    }
}

int daaRecursive(char* s1, char* s2, int len1, int len2) {
    if (len1 == 0 || len2 == 0) {
        printf("Alignment:\n%s\n%s\n", s1, s2);
        return 0;
    }
    
    if (len1 == 1 && len2 == 1) {
        int score = substitutionScore(s1[0], s2[0]);
        printf("Alignment:\n%c\n%c\n", s1[0], s2[0]);
        return score;
    }
    
    int split1 = len1 / 2;
    int split2 = len2 / 2;
    
    int** firstPart = buildDPTable(s1, s2, split1, split2);
    printf("Upper-Left DP Matrix:\n");
    printMatrix(firstPart, split1 + 1, split2 + 1);
    printf("Path:\n");
    backtrackAlignment(s1, s2, split1, split2, firstPart);
    
    int** secondPart = buildDPTable(s1 + split1, s2 + split2, 
                                    len1 - split1, len2 - split2);
    printf("\nLower-Right DP Matrix:\n");
    printMatrix(secondPart, len1 - split1 + 1, len2 - split2 + 1);
    printf("Path:\n");
    backtrackAlignment(s1 + split1, s2 + split2, 
                       len1 - split1, len2 - split2, secondPart);
    
    int partA = firstPart[split1][split2];
    int partB = secondPart[len1 - split1][len2 - split2];
    
    for (int i = 0; i <= split1; i++)
        free(firstPart[i]);
    free(firstPart);
    
    for (int i = 0; i <= len1 - split1; i++)
        free(secondPart[i]);
    free(secondPart);
    
    return partA + partB;
}

int main() {
    char seq1[] = "AGTGATG";
    char seq2[] = "GTTAG";
    int len1 = strlen(seq1);
    int len2 = strlen(seq2);
    
    int result = daaRecursive(seq1, seq2, len1, len2);
    printf("\nMax Score: %d\n", result);
    
    return 0;
}