Python计算FASTA序列中目标20-mers累计计数问题求助

问题背景 #

你手头有两个文件需要处理：

• 常规FASTA文件 single_mapped.fa，其中每条reads长度为120 nt
• 制表符分隔的文本文件 20frequent_20mers.txt，包含10000条20-mer序列及其对应计数，示例内容如下：

AAAAAGTATAGGAGATAGAA 35
AAAAATAGGAGGACTATTCA 26
AAAAATAGGAGGACTATTTA 24
AAAAATAGGAGGCCTATTCA 62

需求说明 #

遍历single_mapped.fa中的每条reads，计算该reads中所有出现在20frequent_20mers.txt里的20-mers的计数累加值。比如某reads包含序列AAAAAGTATAGGAGATAGAA和AAAAATAGGAGGACTATTCA，累加值应为61（35+26）。

现有代码 #

你编写的代码如下：

import csv
from Bio import SeqIO

file2 = open('20frequent_20mers.txt','r')
kmer_list = csv.reader(file2, delimiter='\t')
for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    print(seq_record.id)
    score_fre = 0
    sequence_string = str(seq_record.seq)
    for i in range(0,101):
        seq = sequence_string[i:i+20]
        for row in kmer_list:
            if row[0] == seq:
                score_fre = score_fre + int(row[1])
    print(score_fre)

问题排查 #

你的代码单独运行各循环正常，但组合后失效，核心问题在于**csv.reader返回的是一个迭代器**：

• 第一次遍历kmer_list（也就是第一条reads的循环里），迭代器就已经耗尽了，后续的reads再去遍历kmer_list时，不会有任何数据返回，自然累加值始终为0或者不符合预期。
• 另外，每次检查一个20-mer都要遍历10000条数据，时间复杂度是O(N*M)（N是reads数量，M是每条reads的20-mer数量，这里是101），效率极低，数据量大的时候会非常慢。

高效解决方案 #

我们可以把20-mer和对应的计数预先存入一个字典（查找时间复杂度O(1)），这样既解决了迭代器耗尽的问题，又大幅提升了运行效率：

import csv
from Bio import SeqIO

# 第一步：把高频20-mer加载到字典中，key是序列，value是计数
kmer_counts = {}
with open('20frequent_20mers.txt', 'r') as file2:
    reader = csv.reader(file2, delimiter='\t')
    for seq, count in reader:
        kmer_counts[seq] = int(count)

# 第二步：遍历FASTA序列计算累加值
for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    print(seq_record.id)
    sequence = str(seq_record.seq)
    # 计算所有可能的20-mer，累加对应计数
    total_score = 0
    # 通用写法：长度为L的序列，20-mer的数量是 L - 20 + 1
    for i in range(len(sequence) - 20 + 1):
        current_kmer = sequence[i:i+20]
        # 如果当前20-mer在字典中，就累加计数
        total_score += kmer_counts.get(current_kmer, 0)
    print(total_score)

进一步优化（更简洁的写法） #

可以用生成器表达式结合sum()函数简化累加部分，代码更紧凑：

import csv
from Bio import SeqIO

with open('20frequent_20mers.txt', 'r') as file2:
    kmer_counts = {seq: int(count) for seq, count in csv.reader(file2, delimiter='\t')}

for seq_record in SeqIO.parse("single_mapped.fa", "fasta"):
    sequence = str(seq_record.seq)
    total_score = sum(
        kmer_counts.get(sequence[i:i+20], 0)
        for i in range(len(sequence) - 20 + 1)
    )
    print(f"{seq_record.id}\t{total_score}")

额外说明 #

• 使用with语句打开文件可以自动管理文件资源，避免忘记关闭文件的问题。
• 字典的get()方法可以在20-mer不存在时返回默认值0，不需要额外的判断逻辑。
• 通用的len(sequence)-20+1写法，即使reads长度发生变化也能正常工作，不需要硬编码101。

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