如何用Python从DNA序列6种读码框中识别最长ORF

完整实现代码 #

from Bio import Seq
import regex as re

# 反向互补碱基映射表
ReverseCompDNA = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'}
def reverse_complement(seq):
    """
    输入DNA序列，返回其反向互补序列
    """
    return ''.join([ReverseCompDNA[nuc] for nuc in seq])[::-1]

def find_longest_orf(input_seq):
    startP = re.compile('ATG')
    # 清理输入序列的换行符
    clean_seq = input_seq.replace('\n','')
    # 生成反向互补序列
    rev_comp_seq = reverse_complement(clean_seq)
    # 初始化最长ORF记录：(长度, 起始坐标, 终止坐标, 链方向, 蛋白序列, ORF核酸序列)
    longest_orf = (0,)
    
    # 遍历正向、反向两条链，分别处理所有ORF
    for strand, current_seq in [('+', clean_seq), ('-', rev_comp_seq)]:
        # 匹配所有重叠的ATG起始位点，自动覆盖3种读码框偏移
        for match in startP.finditer(current_seq, overlapped=True):
            start_pos = match.start()
            # 翻译到终止密码子为止
            protein = Seq.Seq(current_seq)[start_pos:].translate(to_stop=True)
            protein_len = len(protein)
            # 仅保留更长的ORF
            if protein_len > longest_orf[0]:
                orf_length = protein_len * 3 + 3
                orf_nuc = current_seq[start_pos: start_pos + orf_length]
                # 反向链的坐标转换为原始正向链的坐标
                if strand == '+':
                    ref_start = start_pos
                    ref_end = start_pos + orf_length
                else:
                    total_len = len(clean_seq)
                    ref_start = total_len - (start_pos + orf_length)
                    ref_end = total_len - start_pos
                # 更新最长ORF记录
                longest_orf = (
                    protein_len,
                    ref_start,
                    ref_end,
                    strand,
                    str(protein),
                    orf_nuc
                )
    return longest_orf

使用示例 #

# 测试用输入序列
input_dna = """ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG
TTACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC"""

result = find_longest_orf(input_dna)
print(f"最长ORF编码氨基酸长度：{result[0]}")
print(f"原始序列起始坐标：{result[1]}")
print(f"原始序列终止坐标：{result[2]}")
print(f"所在链方向：{result[3]}")
print(f"编码蛋白序列：{result[4]}")
print(f"ORF核酸序列：{result[5]}")

关键说明 #

• 自动覆盖6种读码框：正向链3种偏移、反向互补链3种偏移，通过匹配所有ATG起始位点自然实现不同读码框的遍历
• 反向链坐标自动转换为原始输入序列的正向坐标，无需手动换算
• 输出包含完整的ORF信息，可直接用于后续分析

注意：如果输入序列包含N等简并碱基，只需在ReverseCompDNA字典中补充对应的碱基映射规则即可。

内容的提问来源于stack exchange，提问作者Become a Briniac with HR