Python读取FASTA文件序列ID匹配失败问题求助

问题场景 #

处理包含多个条目的FASTA文件时，运行给定的Python代码总会触发"NOTE"提示，无法正常解析rRNA亚基信息。已知目标物种名为Octopus vulgaris，不想修改核心代码，仅通过调整FASTA文件头适配现有逻辑。

FASTA文件示例 #

>AF528955.1 Octopus vulgaris isolate 2 16S ribosomal RNA gene, partial sequence; mitochondrial gene for mitochondrial product
UAACUUAUCUUCUAAGCAAAAAACUUGGUUUAUUCUUCAACUAAACUCAAAAUUAAGGAAGUUAAUACAA
UACUUUAAAUUAAUUUUAUUCCUUGAUCACCCCAACCAAAGUUAUUUACAAAUAAAUUAUAUAUAUACAU
AUAUAUCUAUAUAAUAACA
>AF528954.1 Octopus vulgaris isolate 1 16S ribosomal RNA gene, partial sequence; mitochondrial gene for mitochondrial product
UAACUUAUCUUCUAAGCAAAAAACUUGGUUUAUUCUUCAACUAAACUCAAAAUUAAGGAAGUUAAUACAA
UACUUUAAAUUAAUUUUAUUCCUUGAUCACCCCAACCAAAGUUAUUUACAAAUAAAUUAUAUAUAUACAU
AUAUAUCUAUAUAAUAACA

触发问题的核心代码片段 #

def count_rrna_subunits(full_species_name, seq_filename):
    counts = {}
    with MultiSequenceFile(seq_filename, ignore_missing_file=True) as input_seqs:
        for seq_id, seq_desc, seq in input_seqs:
            # 正则要求seq_id格式为：非数字字符串_数字开头的字符串
            m = re.search(rb'^([^0-9]+)_([0-9][^_]+)', seq_id)
            if m:
                species = m.group(1).replace(b'_', b' ')
                subunit = m.group(2)
                if species == full_species_name.encode('ascii'):
                    counts[subunit] = counts.get(subunit, 0) + 1
            else:
                print("NOTE: could not parse rRNA subunits information for entry:\n{}Reference rRNA subunit listing may be unavailable in miRTrace reports generated using this database. Actual read mapping not affected.\n".format(seq_id), file=sys.stderr)
    return counts

问题根源 #

代码中的正则表达式rb'^([^0-9]+)_([0-9][^_]+)'对seq_id（FASTA头中>后第一个空格前的内容）有严格格式要求：

1. 必须包含下划线_
2. 下划线前是纯非数字字符（用来提取物种名，替换下划线为空格后要匹配Octopus vulgaris）
3. 下划线后以数字开头（用来提取rRNA亚基，比如16S）

当前FASTA的seq_id是AF528955.1这类编号，既没有下划线，也不符合上述格式，导致正则匹配失败，触发NOTE提示。

解决方案：修改FASTA文件头格式 #

单条修改示例 #

将原FASTA头：

>AF528955.1 Octopus vulgaris isolate 2 16S ribosomal RNA gene...

修改为：

>Octopus_vulgaris_16S AF528955.1 Octopus vulgaris isolate 2 16S ribosomal RNA gene...

新的seq_idOctopus_vulgaris_16S完全符合正则要求：

• 下划线前的Octopus_vulgaris替换下划线为空格后，正好匹配目标物种名Octopus vulgaris
• 下划线后的16S以数字开头，会被识别为rRNA亚基

批量修改命令（Linux/macOS） #

如果FASTA文件条目较多，用sed命令批量处理：

sed -E 's/^>([A-Z0-9]+\.[0-9]) (Octopus vulgaris) isolate [0-9]+ ([0-9]+S)/>\2_\3 \1 \2 isolate \3/' input.fasta > output.fasta

该命令会自动将所有符合格式的条目头替换为符合要求的结构，生成新的FASTA文件output.fasta。

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