Python解析GBFF文件并提取CDS序列至FASTA的方法

GBFF文件的结构组成

GBFF是GenBank数据库的基础交换格式，每个独立条目以双斜杠 // 分隔。核心区块包含：

• LOCUS：记录序列长度、分子类型（DNA/RNA）、物种及提交日期。
• DEFINITION：提供该序列的生物学功能摘要，常作为唯一标识键。
• FEATURES：描述基因产物与调控元件，其中CDS（Coding Sequence）标注翻译为蛋白质的开放阅读框，并附带起始/终止坐标。
• ORIGIN：实际碱基排列区，通常采用每10个字符添加行号索引的排版方式，需清洗后才能转为连续字符串。

使用Python解析GBFF提取CDS的步骤

第1步：初始化变量并读取DEFINITION字段

目的：建立当前条目的标识映射，确保后续坐标仅作用于所属序列。
入口：迭代器逐行读取GBFF文本流
动作：匹配以 DEFINITION 开头的行，剔除尾部句点与换行符后赋值给临时字典键
重要参数：文件编码需设为 utf-8，防止含特殊符号的物种名乱码
完成状态：current_def 成功捕获标题，准备好接收后续坐标

第2步：捕获FEATURES中的CDS位置信息

目的：获取待提取片段的数值边界，跳过蛋白质、重复序列等其他注释
入口：匹配 CDS 关键字所在行
动作：使用正则表达式提取行内所有数字，将首尾组合为闭区间元组存入列表
重要参数：需忽略 complement() 或 join() 包裹的非线性坐标（基础脚本暂作简化处理）
完成状态：current_cds_loc 列表中准确记录 (start, end) 对

第3步：遍历ORIGIN区块提取纯碱基序列

目的：将格式化换行的碱基重组为可直接切片的一维字符串
入口：检测到 ORIGIN 头部标志进入读块状态
动作：持续收集直至遇见结束标记 //，期间剥离行首数字索引、制表符与多余空格，统一转为大写
重要参数：re.sub(r'\d+\s*', '', line) 用于安全去噪
完成状态：full_seq 拼接完毕，长度与LOCUS声明一致

第4步：按坐标截取CDS并输出FASTA文件

目的：根据特征表标注的物理位置提取有效转录单元
入口：ORIGIN区块闭合瞬间触发写入逻辑
动作：执行 full_seq[start-1 : end] 完成闭开区间转换，合并多段提取结果后，按 >Title\nSequence 模板追加至目标文件
重要参数：输出文件需预留足够磁盘空间，建议按批次分卷防内存峰值
完成状态：终端打印提取计数，输出文件夹生成合法FASTA后缀文件

import re

INPUT_GB = "./input.gbff"
OUTPUT_FASTA = "./output.fasta"

def parse_gbff_to_fasta(input_path, output_path):
    with open(input_path, 'r', encoding='utf-8') as f_in, \
         open(output_path, 'w', encoding='utf-8') as f_out:

        current_def = None
        current_cds_loc = []
        origin_lines = []
        in_origin = False

        for line in f_in:
            line = line.strip()
            if not line:
                continue

            # 识别定义行，提取标题作为键
            if line.startswith('DEFINITION'):
                title = line.replace('DEFINITION ', '').replace('.', '').strip()
                current_def = title
                current_cds_loc = []
                origin_lines = []
                in_origin = False
                print(f"解析条目: {current_def}")

            # 捕获 CDS 特性及其位置
            elif line.startswith('CDS'):
                loc_str = line.split('CDS', 1)[1].strip()
                coords = re.findall(r'(\d+)', loc_str)
                if len(coords) >= 2:
                    current_cds_loc.append((int(coords[0]), int(coords[-1])))

            # 进入 ORIGIN 区域，拼接纯碱基序列
            elif line.startswith('ORIGIN'):
                in_origin = True
                origin_lines = []
            elif in_origin:
                if line == '//':
                    in_origin = False
                    # 清理空位、制表符和行首数字索引，合并为完整序列
                    seq_parts = []
                    for l in origin_lines:
                        cleaned = re.sub(r'\d+\s*', '', l.replace('\t', ''))
                        if cleaned:
                            seq_parts.append(cleaned.upper())
                    full_seq = ''.join(seq_parts)

                    # 按之前记录的坐标截取 CDS
                    if current_def and current_cds_loc:
                        cds_seq = []
                        for start, end in current_cds_loc:
                            # GBFF坐标为闭区间，Python切片为左闭右开，故减1
                            seg = full_seq[start - 1 : end]
                            cds_seq.append(seg)
                        
                        if cds_seq:
                            cds_combined = ''.join(cds_seq)
                            f_out.write(f">{current_def}\n{cds_combined}\n")
                else:
                    origin_lines.append(line)

if __name__ == "__main__":
    parse_gbff_to_fasta(INPUT_GB, OUTPUT_FASTA)
    print("提取完成，已保存至指定 FASTA 文件。")

注意事项与排错建议

在实际生物信息学流水线中，需注意以下细节以避免结果偏差：

• 方向性与反向互补：部分CDS位于负链且标记为 complement()。本脚本按正向坐标直接截取，若需严格还原氨基酸编码链，应在提取后调用序列反转补工具处理。
• 复杂断点连接：真核生物基因常存在内含子，坐标形式为 join(100..200, 500..600)。基础正则会误判首尾数字，建议改用成熟的 Bio.SeqRecord 对象解析 features 属性。
• 乱码与编码：早期GenBank提交可能混用ISO-8859-1编码，读取失败时请尝试切换 encoding='latin1' 或先运行一次全局转码。

以上就是Python解析GBFF文件并提取CDS序列的详细内容，更多关于Python生物信息处理与基因数据格式转换的资料请关注本站其它相关文章！