基于GFF3与FASTA提取对应区间序列并合并至Pandas DataFrame

现有数据 #

GFF3导入后得到的DataFrame #

seq_id  source  type     start       end score strand phase      attributes
1      ctg.s1.000000F_arrow  prokka  gene      56.0     244.0     .      +     .  NHDIEHOJ_00001
3      ctg.s1.000000F_arrow  prokka  gene     902.0    2167.0     .      -     .  NHDIEHOJ_00002
5      ctg.s1.000001F_arrow  prokka  gene    2363.0    2905.0     .      -     .  NHDIEHOJ_00003
7      ctg.s1.000003F_arrow  prokka  gene    2916.0    3947.0     .      -     .  NHDIEHOJ_00004
9      ctg.s2.000000F_arrow  prokka  gene    4353.0    5174.0     .      +     .  NHDIEHOJ_00005

seq_id列的唯一值 #

ctg.s1.000000F_arrow
ctg.s1.000001F_arrow
ctg.s1.000003F_arrow
ctg.s2.000000F_arrow

已转换为字典的FASTA数据 #

FASTA文件已转换为Python字典，键为FASTA头部的标识（即上述seq_id），值为对应的大写序列：

>ctg.s1.000000F_arrow
CCGGAACATCGTCCTCATCCGCAAAGTCGAGCTCTGCCTCGATCATTGCACGCGAATGGGTCAGCCGTCGGGCCCAACCG
GCATAGAGTGCGGACTGTCCGCCACCGGACTGCTCTATGGCGAGACGACGCTGCATTTCCGTTTCTGCCGCAATCAGGTC
>ctg.s1.000001F_arrow
ACGCCGGCTGCAACTTTGAGAAGATGTGGCGATGTCGACCGCTGCATCCCGCCCTTCTCTGCAGAATTTTCCAGCTGTCC
GAGGACATTGGCAAAAAGGCCCTTGGAATCCTTGCGGCTCATTCTCTCCCCCATGCCTTCCAGAAGAGGCCCTCGAGTTC
>ctg.s1.000003F_arrow
GGCGCTGGTTTTCCCCGACACCTCGCCGCGCGGCGAGGGCGTGGCTGACGACGAGGCTTATGATCTCGGTCAGGGTGCGG
GCTTCTATGTCAATGCGACGCAGAAGCCCTGGTCGCCGCACTATCGCATGTATGATTATATCGTCACCGAATTGCCCGCC
>ctg.s2.000000F_arrow
GCGCTCGACGGCATGCCCGTACGCGGCCGATCCTGCGCCGCTTCCTTAACCTTAGCTGCGGATGGAAAGTCGTCCTCGGA
GTTCGGCTCGCAAACGCTTTCGAGCGCGCAATTGACGACGATGTCGTACCCAACTTAGATCGCCGAACGCCATGAGGTCG

需求 #

为上述DataFrame新增sequence列：

• 每一行根据seq_id匹配FASTA字典中的对应序列
• 按照start和end提取区间序列（注意GFF3坐标为1-based，需转换为Python的0-based索引）
• 若strand为-，需对提取的序列做反向互补处理
• sequence列长度需与end - start + 1（区间实际长度）一致

最终目标DataFrame示例：

seq_id  source  type     start       end score strand phase      attributes    sequence
1      ctg.s1.000000F_arrow  prokka  gene      56.0     244.0     .      +     .  NHDIEHOJ_00001    CCGGAACATCGTCCTCATCCG...
3      ctg.s1.000000F_arrow  prokka  gene     902.0    2167.0     .      -     .  NHDIEHOJ_00002    CAAGGACATCGTGATCAATTC...
5      ctg.s1.000001F_arrow  prokka  gene    2363.0    2905.0     .      -     .  NHDIEHOJ_00003    TCGCCGCGCGGCGAGTGATTA...
7      ctg.s1.000003F_arrow  prokka  gene    2916.0    3947.0     .      -     .  NHDIEHOJ_00004    TCGAGCGCGCAATTGACGACG...
9      ctg.s2.000000F_arrow  prokka  gene    4353.0    5174.0     .      +     .  NHDIEHOJ_00005    AGATCGCCGAACGCCATATTT...

解决方案代码 #

import pandas as pd

# 假设已导入的GFF DataFrame命名为gff_df
# 假设FASTA字典命名为fasta_dict

# 定义反向互补函数
def reverse_complement(seq):
    complement_map = {'A':'T', 'T':'A', 'C':'G', 'G':'C'}
    return ''.join([complement_map[base] for base in reversed(seq)])

# 定义提取序列的行处理函数
def extract_sequence(row):
    seq_id = row['seq_id']
    # 转换GFF3的1-based坐标为Python字符串的0-based索引
    start_idx = int(row['start']) - 1
    end_idx = int(row['end'])
    strand = row['strand']
    
    # 从FASTA字典获取完整序列
    full_seq = fasta_dict[seq_id]
    # 提取目标区间序列
    region_seq = full_seq[start_idx:end_idx]
    
    # 负链序列做反向互补处理
    if strand == '-':
        region_seq = reverse_complement(region_seq)
    
    return region_seq

# 为DataFrame添加sequence列
gff_df['sequence'] = gff_df.apply(extract_sequence, axis=1)

# 查看结果
print(gff_df)

关键说明 #

1. 坐标转换：GFF3采用1-based坐标，Python字符串为0-based，因此start需减1，end直接使用（切片[start:end]包含start到end-1，刚好对应1-based的start到end区间）
2. 负链处理：反向互补是基因组分析中负链序列的标准处理逻辑，确保序列方向与基因转录方向一致
3. 长度匹配：提取后的序列长度为end - start（0-based切片），对应1-based的end - start + 1，完全匹配区间实际长度

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