如何优化Slurm上运行的Python FASTA解析脚本提升运行速度

你这个问题属于Python代码性能优化+生物信息学脚本开发范畴，在对应技术板块发帖是完全匹配的。

• 字符串拼接的时间复杂度爆炸：Python中字符串是不可变类型，你逐行读取序列时用dictionnary[seqid] += record做拼接，每次执行都会重新分配内存、复制该条序列之前已经读入的全部内容，单条序列越长，拼接的耗时会呈平方级增长。小型测试基因组序列短感知不明显，碰到哺乳动物级别几百Mb甚至上Gb的染色体序列，这个写法的耗时会直接涨到数小时到数天级别。
• 无意义的控制台IO拖慢速度：解析序列ID时你对每个seqid都执行print(seqid)，后续还直接print(all_len)打印全量序列长度列表。哺乳动物基因组的contig/scaffold数量通常在数万到数十万级别，高频的控制台输出本身就会占用大量运行资源。
• 冗余内存占用与重复计算：你把所有完整基因组序列全量加载到内存中存储，仅为了后续统计序列长度和N碱基数量，哺乳动物全基因组数据加载到内存会占用数GB到十几GB空间，触发内存交换时会进一步拖慢速度；而且存完全量序列后再遍历算长度、算N碱基数量，属于完全没必要的重复操作，逐行读取时就可以同步完成统计。

优化逻辑：逐行读取时直接统计序列长度、N碱基数量，不做逐行字符串拼接，去掉无意义的全量打印，内存占用降低90%以上，处理哺乳动物基因组的耗时可以从数天压缩到数秒到数十秒级别，同时兼容无>开头标识的序列输入场景：

from collections import OrderedDict
import argparse

parser = argparse.ArgumentParser(description = "N50 parser")
parser.add_argument("-i", "--input", action="store", dest="input", required=True,
                    help="Input file with sequences")
parser.add_argument("-o", "--output", action="store", dest="output",
                    help="output file for N50/L50 result")
parser.add_argument("-o2", "--output2", action="store", dest="output2",
                    help="output file for sequence length list")
args = parser.parse_args()


def read_fasta(fasta_file):
    seq_length_dict = OrderedDict()
    description_dict = OrderedDict()
    total_N = 0
    current_seqid = None
    # 兼容没有>开头标识的FASTA场景，初始化默认ID
    default_id = "seq_0"
    seq_length_dict[default_id] = 0
    description_dict[default_id] = ""
    current_seqid = default_id

    with open(fasta_file, 'r') as infile:
        for line in infile:
            record = line.strip()
            if not record:
                continue
            if record[0] == '>':
                seqid = record.split(" ")[0][1:]
                seq_length_dict[seqid] = 0
                toto = record.split(" ", 1)
                description_dict[seqid] = toto[1] if len(toto)>=2 else ""
                current_seqid = seqid
                continue
            # 逐行直接累加长度、统计N数量，不存全序列
            seq_line_len = len(record)
            seq_length_dict[current_seqid] += seq_line_len
            total_N += record.count("N")
    # 如果文件里存在>开头的行，删掉初始化的空默认ID
    if len(seq_length_dict) > 1 and seq_length_dict[default_id] == 0:
        del seq_length_dict[default_id]
        del description_dict[default_id]
    return seq_length_dict, description_dict, total_N

seq_length_dict, descriptdict, N_number = read_fasta(args.input)
total_length = sum(seq_length_dict.values())
all_len = sorted(seq_length_dict.values(), reverse=True)

L50 = 0
N50 = 0
if total_length > 0:
    acum = 0
    half_total = total_length / 2
    for idx, contig_len in enumerate(all_len):
        acum += contig_len
        if acum >= half_total:
            L50 = idx + 1
            N50 = contig_len
            break

print(f"Total N count: {N_number}")
print(f"Total assembly length: {total_length}")
print(f"The L50 is {L50}")
print(f"The N50 is {N50}")

with open(args.output, 'w') as outfile:
    outfile.write(f"L50\t{L50}\n")
    outfile.write(f"N50\t{N50}\n")

with open(args.output2, "w") as file:
    for key, value in seq_length_dict.items():
        file.write(f"{key} : {value}\n")

如果后续你确实需要提取完整序列做其他分析，不要用+=逐行拼接字符串，把每条序列的行存在列表里，最后用''.join(seq_line_list)一次性合并成完整字符串，速度会比逐行拼接快几个数量级。

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