如何用Python高效统计4G级FASTA文件的字符频率？

你当前的代码耗时较长，核心原因是逐行遍历+逐个字符的Python循环本身开销大，再加上频繁的进度打印、字符串处理操作，进一步拖慢了速度。以下是几种更高效的解决方案，按速度优先级排序：

1. 用命令行工具awk（最快，原生C实现） #

直接用awk处理，无需编写Python代码，速度比Python快数倍甚至一个数量级：

time awk '!/^>/{gsub(/[^AGTCN]/,""); for(i=1;i<=length($0);i++) c[toupper(substr($0,i,1))]++} END{for(k in c) print k,c[k]}' dna_sequences.fa

原理：跳过所有>开头的行，过滤掉非目标碱基的字符，逐个统计大写后的碱基出现次数，最后输出结果。

2. 用BioPython专业库（生物信息学场景最优） #

BioPython的SeqIO模块专门针对FASTA等生物序列文件做了优化，自动处理注释行，代码简洁且效率高：

import time
from Bio import SeqIO

def count_bases(file_path):
    counts = {"A": 0, "G": 0, "T": 0, "C": 0, "N": 0}
    begin = time.time()

    for record in SeqIO.parse(file_path, "fasta"):
        seq = str(record.seq).upper()
        counts['A'] += seq.count('A')
        counts['G'] += seq.count('G')
        counts['T'] += seq.count('T')
        counts['C'] += seq.count('C')
        counts['N'] += seq.count('N')
    
    end = time.time()
    print(f"总耗时: {end - begin:.2f}s")
    return counts

if __name__ == "__main__":
    result = count_bases("dna_sequences.fa")
    print(result)

安装BioPython：pip install biopython

3. Python批量字符串统计（内存足够时首选） #

利用Python字符串的count方法（底层C实现），一次性读取内容后过滤注释行，批量统计：

import time
import re

def count_bases(file_path):
    counts = {}
    target_chars = ['A', 'G', 'T', 'C', 'N']
    begin = time.time()

    with open(file_path, 'rt') as f:
        # 过滤所有>开头的行
        content = re.sub(r'^>.*\n?', '', f.read(), flags=re.MULTILINE)
        content = content.upper()
        for char in target_chars:
            counts[char] = content.count(char)
    
    end = time.time()
    print(f"总耗时: {end - begin:.2f}s")
    return counts

if __name__ == "__main__":
    result = count_bases("dna_sequences.fa")
    print(result)

注意：此方法需要足够内存（4G文件读取后内存占用约4-6G），内存不足时慎用。

4. 分块读取优化（内存有限场景） #

如果内存不够，采用分块读取的方式，减少单次内存占用，同时用Counter批量更新统计：

import time
from collections import Counter

def count_bases(file_path):
    counts = Counter({"A":0, "G":0, "T":0, "C":0, "N":0})
    block_size = 1024 * 1024 * 64  # 64MB块，可根据内存调整
    skip_next_line = False

    begin = time.time()
    with open(file_path, 'rt') as f:
        while True:
            block = f.read(block_size)
            if not block:
                break
            lines = block.split('\n')
            # 处理上一块遗留的需跳过的行
            if skip_next_line:
                lines.pop(0)
                skip_next_line = False
            # 遍历完整行
            for line in lines[:-1]:
                if line.startswith(">"):
                    skip_next_line = True
                    continue
                # 只统计目标碱基
                counts.update(c.upper() for c in line.strip() if c.upper() in counts)
            # 处理块内最后一行（可能不完整）
            last_line = lines[-1]
            if last_line.startswith(">"):
                skip_next_line = True
            else:
                counts.update(c.upper() for c in last_line.strip() if c.upper() in counts)
    
    end = time.time()
    print(f"总耗时: {end - begin:.2f}s")
    return counts

if __name__ == "__main__":
    result = count_bases("dna_sequences.fa")
    print(result)

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