如何从FASTA文件中筛选可翻译为蛋白质的有效氨基酸序列？

需求明确 #

从FASTA文件里挑出有效序列，删掉无效的，同时给有效序列提取完整读框后存到新文件：

• 无效序列：要么不含起始氨基酸M，要么终止密码子*占比超过你设定的阈值（比如5%）
• 完整读框：从序列里第一个M开始，到第一个*结束（不含*；要是序列里没*，就保留从M开始的全部内容）

实现脚本（Python） #

def parse_fasta(file_path):
    """把FASTA文件拆成序列ID和对应的完整氨基酸序列"""
    sequences = {}
    current_id = None
    current_seq = []
    with open(file_path, 'r') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            if line.startswith('>'):
                if current_id is not None:
                    sequences[current_id] = ''.join(current_seq)
                current_id = line[1:]
                current_seq = []
            else:
                current_seq.append(line)
        # 处理最后一条没写完的序列
        if current_id is not None:
            sequences[current_id] = ''.join(current_seq)
    return sequences

def filter_and_extract(sequences, stop_codon_threshold=0.05):
    """筛选有效序列，同时提取完整读框"""
    valid_sequences = {}
    for seq_id, seq in sequences.items():
        # 先看有没有起始M
        if 'M' not in seq:
            continue
        # 算终止密码子占比
        stop_count = seq.count('*')
        stop_ratio = stop_count / len(seq) if len(seq) > 0 else 0
        if stop_ratio > stop_codon_threshold:
            continue
        # 提取读框：从第一个M到第一个*（不含*）
        start_idx = seq.find('M')
        end_idx = seq.find('*', start_idx)
        if end_idx == -1:
            # 没找到终止密码子，就留从M开始的所有内容
            full_frame = seq[start_idx:]
        else:
            full_frame = seq[start_idx:end_idx]
        # 确保提取的序列不是空的
        if len(full_frame) > 0:
            valid_sequences[seq_id] = full_frame
    return valid_sequences

def write_fasta(sequences, output_path):
    """把处理好的序列写入FASTA，每行最多80个字符，符合标准格式"""
    with open(output_path, 'w') as f:
        for seq_id, seq in sequences.items():
            f.write(f'>{seq_id}\n')
            # 分块写，避免一行太长
            for i in range(0, len(seq), 80):
                f.write(f'{seq[i:i+80]}\n')

if __name__ == '__main__':
    # 替换成你的实际文件路径
    input_file = 'example.fasta'
    output_file = 'valid_sequences.fasta'
    # 终止密码子占比阈值，自己调，比如允许10%就改成0.1
    stop_threshold = 0.05  # 5%
    
    # 跑流程
    fasta_seqs = parse_fasta(input_file)
    valid_seqs = filter_and_extract(fasta_seqs, stop_threshold)
    write_fasta(valid_seqs, output_file)
    
    print(f"搞定！有效序列已经存到 {output_file} 里了")

使用步骤 #

1. 把上面的代码存成filter_fasta.py文件
2. 修改input_file和output_file为你自己的文件路径
3. 根据需求调整stop_threshold（终止密码子占比的上限）
4. 打开终端跑命令：python filter_fasta.py

针对示例文件的处理结果 #

示例里的三条序列：

• abc：有M，没有终止密码子 → 有效，提取从第一个M开始的全部序列
• def：虽然有M，但终止密码子占比超过5% → 被判定为无效，直接删掉
• ghi：有M，终止密码子占比不到5% → 有效，提取从第一个M到第一个*之前的部分

最终生成的valid_sequences.fasta会包含abc和ghi的完整读框序列。

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