如何修改C语言序列比对代码以加载指定文件并实现半全局比对？

需求说明 #

现有一段实现双序列全局比对的C语言动态规划代码，需修改为支持加载两个独立文件（分别对应SARS-CoV-2参考基因组和Nanopore读段），并适配半全局比对逻辑（用于短读段与长参考序列的局部区域比对）。原代码仅支持单输入文件，运行时提示Usage: align < input file >。

关键修改步骤 #

• 调整命令行参数逻辑：支持传入两个文件名，修改参数检查逻辑，更新用法提示为Usage: align <ref_genome_file> <nanopore_read_file>
• 修复序列存储溢出问题：将原代码中固定长度为1000的序列数组，替换为宏定义的MAXLENGTH_A和MAXLENGTH_B，确保能容纳完整的基因组和读段序列
• 实现双文件读取：分别打开两个输入文件，读取各自的序列内容（假设每个文件仅包含一行无注释的碱基序列）
• 适配半全局比对规则：修改动态规划矩阵的初始化逻辑，允许参考序列或读段两端的gap不扣分；调整回溯起点为最后一行/列的最高得分位置，而非全局比对的右下角
• 修复内存分配硬编码：将内存分配中的硬编码数值替换为宏定义，提升代码可维护性
• 修正argv索引错误：原代码中错误使用argv[2]，改为对应正确的参数索引

修改后的完整代码 #

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

#define GAP -2
#define MATCH 5
#define MISMATCH -3

#define MAXLENGTH_A 29904  // 对应SARS-CoV-2参考基因组长度+1
#define MAXLENGTH_B 1247   // 对应Nanopore读段长度+1

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;
}

char Tmax(int A, int B, int C)
{
    if (A > B && A > C) return 'D';
    else if (B > A && B > C) return 'L';
    else return 'U';
}

int m(char p, char q)
{
    if (p == q) return MATCH;
    else return MISMATCH;
}

void append(char *st, int L, char c)
{
    int i;
    for (i = L; i > 0; i--)
        st[i] = st[i-1];
    st[L+1] = '\0';
    st[0] = c;
}

// 读取文件中的序列（假设文件仅含一行碱基序列）
void read_sequence(const char *filename, char *seq, int max_len)
{
    FILE *fp = fopen(filename, "r");
    if (!fp) {
        fprintf(stderr, "Error: Cannot open file %s\n", filename);
        exit(1);
    }
    // 读取整行序列，忽略换行符
    if (!fgets(seq, max_len, fp)) {
        fprintf(stderr, "Error: Failed to read sequence from %s\n", filename);
        fclose(fp);
        exit(1);
    }
    // 移除换行符（如果存在）
    seq[strcspn(seq, "\n")] = '\0';
    fclose(fp);
}

int main(int argc, char **argv)
{
    char A[MAXLENGTH_A];
    char B[MAXLENGTH_B];
    char RA[MAXLENGTH_A + MAXLENGTH_B];  // 足够存储比对后的序列
    char RM[MAXLENGTH_A + MAXLENGTH_B];
    char RB[MAXLENGTH_A + MAXLENGTH_B];
    int N, M, L;
    int i, j;

    // 动态分配DP矩阵
    int **S = (int**)malloc(sizeof(int*) * MAXLENGTH_A);
    for (i = 0; i < MAXLENGTH_A; i++)
        S[i] = (int*)malloc(sizeof(int) * MAXLENGTH_B);
    char **T = (char**)malloc(sizeof(char*) * MAXLENGTH_A);
    for (i = 0; i < MAXLENGTH_A; i++)
        T[i] = (char*)malloc(sizeof(char) * MAXLENGTH_B);

    // 检查命令行参数
    if (argc != 3) {
        fprintf(stderr, "Usage: align <SARS-CoV-2_ref_genome.txt> <Nanopore.txt>\n");
        exit(1);
    }

    // 读取两个文件的序列
    read_sequence(argv[1], A, MAXLENGTH_A);
    read_sequence(argv[2], B, MAXLENGTH_B);

    printf("Sequence A length: %zu\n", strlen(A));
    printf("Sequence B length: %zu\n", strlen(B));

    N = strlen(A);
    M = strlen(B);

    // 半全局比对初始化：允许两端gap不扣分
    for (i = 0; i <= N; i++) {
        S[i][0] = 0;
        T[i][0] = 'U';
    }
    for (j = 0; j <= M; j++) {
        S[0][j] = 0;
        T[0][j] = 'L';
    }
    T[0][0] = 'D';

    // 填充DP矩阵
    for (i = 1; i <= N; i++) {
        for (j = 1; j <= M; j++) {
            int match_score = S[i-1][j-1] + m(A[i-1], B[j-1]);
            int gap_A = S[i][j-1] + GAP;
            int gap_B = S[i-1][j] + GAP;
            S[i][j] = max(match_score, gap_A, gap_B);
            T[i][j] = Tmax(match_score, gap_A, gap_B);
        }
    }

    // 找到半全局比对的最高得分位置（最后一行或最后一列）
    int max_score = S[N][0];
    int end_i = N, end_j = 0;
    // 检查最后一行
    for (j = 1; j <= M; j++) {
        if (S[N][j] > max_score) {
            max_score = S[N][j];
            end_i = N;
            end_j = j;
        }
    }
    // 检查最后一列
    for (i = 1; i <= N; i++) {
        if (S[i][M] > max_score) {
            max_score = S[i][M];
            end_i = i;
            end_j = M;
        }
    }

    printf("The score of the semi-global alignment is: %d\n", max_score);

    // 回溯生成比对结果
    i = end_i;
    j = end_j;
    L = 0;
    RA[0] = '\0';
    RB[0] = '\0';
    RM[0] = '\0';

    while (i != 0 || j != 0) {
        if (T[i][j] == 'D') {
            append(RA, L, A[i-1]);
            append(RB, L, B[j-1]);
            append(RM, L, (A[i-1] == B[j-1]) ? '|' : '*');
            i--;
            j--;
        } else if (T[i][j] == 'L') {
            append(RA, L, '-');
            append(RB, L, B[j-1]);
            append(RM, L, ' ');
            j--;
        } else if (T[i][j] == 'U') {
            append(RA, L, A[i-1]);
            append(RB, L, '-');
            append(RM, L, ' ');
            i--;
        }
        L++;
    }

    printf("\nAlignment result:\n");
    printf("%s\n", RA);
    printf("%s\n", RM);
    printf("%s\n", RB);

    // 释放动态分配的内存
    for (i = 0; i < MAXLENGTH_A; i++) {
        free(S[i]);
        free(T[i]);
    }
    free(S);
    free(T);

    return 0;
}

使用说明 #

编译代码：

gcc -o align align.c

运行程序：

./align SARS-CoV-2_ref_genome.txt Nanopore.txt

内容的提问来源于stack exchange，提问作者McKale Grant