SARS-CoV-2参考基因组与Nanopore读段半全局比对代码报错求助

问题背景 #

用以下动态规划代码对29903bp的SARS-CoV-2参考基因组和1246bp的Nanopore样本首条读段做半全局比对时，终端输出报错：Usage: align <input file>。

原代码 #

import sys
import numpy as np
GAP = -2
MATCH = 5
MISMATCH = -3
MAXLENGTH_A = 29904
MAXLENGTH_B = 1247

# insert sequence files
A = open("SARS-CoV-2 reference genome.txt", "r")
B = open("Nanopore.txt", "r")

def max(A, B, C):
    if (A >= B and A >= C):
        return A
    elif (B >= A and B >= C):
        return B
    else:
        return C
def Tmax(A, B, C):
    if (A > B and A > C):
        return 'D'
    elif (B > A and B > C):
        return 'L'
    else:
        return 'U'
def m(p, q):
    if (p == q):
        return MATCH
    else:
        return MISMATCH
def append(st, c):
    return c + "".join(i for i in st)
if __name__ == "__main__":
    if (len(sys.argv) != 2):
        print("Usage: align <input file>")
        sys.exit()
    if (not os.path.isfile(sys.argv[1])):
        print("input file not found.")
        sys.exit()
    S = np.empty([MAXLENGTH_A, MAXLENGTH_B], dtype = int)
    T = np.empty([MAXLENGTH_A, MAXLENGTH_B], dtype = str)
    with open(sys.argv[1], "r") as file:
        A = str(A.readline())[:-1]
        B = str(B.readline())[:-1]
        print("Sequence A:",A)
        print("Sequence B:",B)
        N = len(A)
        M = len(B)
        S[0][0] = 0
        T[0][0] = 'D'
        for i in range(0, N + 1):
            S[i][0] = GAP * i
            T[i][0] = 'U'
        for i in range(0, M + 1):
            S[0][i] = GAP * i
            T[0][i] = 'L'
        for i in range(1, N + 1):
            for j in range(1, M + 1):
                S[i][j] = max(S[i-1][j-1]+m(A[i-1],B[j-1]),S[i][j-1]+GAP,S[i-1][j]+GAP)
                T[i][j] = Tmax(S[i-1][j-1]+m(A[i-1],B[j-1]),S[i][j-1]+GAP,S[i-1][j]+GAP)
        print("The score of the alignment is :",S[N][M])
        i, j = N, M
        RA = RB = RM = ""
        while (i != 0 or j != 0):
            if (T[i][j]=='D'):
                RA = append(RA,A[i-1])
                RB = append(RB,B[j-1])
                if (A[i-1] == B[j-1]):
                    RM = append(RM,'|')
                else:
                    RM = append(RM,'*')
                i -= 1
                j -= 1
            elif (T[i][j]=='L'):
                RA = append(RA,'-')
                RB = append(RB,B[j-1])
                RM = append(RM,' ')
                j -= 1
            elif (T[i][j]=='U'):
                RA = append(RA,A[i-1])
                RB = append(RB,'-')
                RM = append(RM,' ')
                i -= 1
        print(RA)
        print(RM)
        print(RB)

报错原因分析 #

1. 命令行参数逻辑矛盾：代码强制要求用户传入一个输入参数，但实际逻辑是直接读取固定文件名的两个序列文件，完全没用到该参数，属于冗余检查导致的错误。
2. 未导入os模块：代码中使用了os.path.isfile()函数，但开头未导入os模块，即使过了参数检查，后续也会触发导入错误。
3. 文件读取逻辑错误：
   • 变量名冲突：先将A、B定义为文件对象，后续又覆盖为字符串，导致读取逻辑混乱；
   • 仅读取第一行：用readline()只读取文件第一行，而基因组FASTA文件通常包含多行序列，会丢失大部分碱基信息。
4. 全局比对而非半全局：代码当前实现的是全局比对（边界初始化时给gap计分），但需求是半全局比对，边界条件不符合要求。

修复步骤 #

1. 移除冗余的命令行参数检查 #

删掉以下代码块，因为代码不需要用户传入额外输入文件：

if (len(sys.argv) != 2):
    print("Usage: align <input file>")
    sys.exit()
if (not os.path.isfile(sys.argv[1])):
    print("input file not found.")
    sys.exit()
with open(sys.argv[1], "r") as file:

2. 导入缺失的os模块 #

在开头补充导入：

import sys
import os
import numpy as np

3. 修正序列读取逻辑 #

添加FASTA文件读取函数，正确拼接多行序列：

def read_fasta(filename):
    with open(filename, 'r') as f:
        seq = ''
        for line in f:
            if not line.startswith('>'):
                seq += line.strip()
    return seq

# 读取序列（替换原来的文件打开代码）
seq_ref = read_fasta("SARS-CoV-2 reference genome.txt")
seq_read = read_fasta("Nanopore.txt")

后续代码中把所有A、B替换为seq_ref、seq_read，避免变量名冲突。

4. 修改为半全局比对的边界条件 #

以允许短读段（seq_read）两端自由gap为例，调整初始化和回溯逻辑：

N = len(seq_ref)
M = len(seq_read)

# 半全局比对初始化：短序列两端gap不计分
S[0][0] = 0
T[0][0] = 'D'
for i in range(0, N + 1):
    S[i][0] = 0  # 短序列开头gap不计分
    T[i][0] = 'U'
for j in range(0, M + 1):
    S[0][j] = 0  # 短序列结尾gap不计分
    T[0][j] = 'L'

# 动态规划计算分数（替换A、B为seq_ref、seq_read）
for i in range(1, N + 1):
    for j in range(1, M + 1):
        match_score = S[i-1][j-1] + m(seq_ref[i-1], seq_read[j-1])
        gap_left = S[i][j-1] + GAP
        gap_up = S[i-1][j] + GAP
        S[i][j] = max(match_score, gap_left, gap_up)
        T[i][j] = Tmax(match_score, gap_left, gap_up)

# 半全局比对找最大分数的回溯起点（最后一行的最大值位置）
max_score = np.max(S[N, :])
j = np.argmax(S[N, :])
i = N
print("The score of the alignment is :", max_score)

5. 其他优化 #

• 避免覆盖内置函数：自定义的max函数覆盖了Python内置的max，建议重命名为max_three；
• 内存优化：对于大序列，可使用滚动数组减少内存占用（仅保留当前行和上一行）。

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