从 fasta 文件估计 Biopython 中的字母表答案

【问题标题】：Estimate Alphabet in Biopython from fasta file从 fasta 文件估计 Biopython 中的字母表
【发布时间】：2017-01-11 10:45:00
【问题描述】：

我正在寻找一种方法来读取 Biopython 中的 .fasta 文件，并估计我们是否正在处理 DNA、RNA 或蛋白质。到目前为止，我读到的数据是这样的：

with open('file.fasta', 'r') as f:
    for seq in sio.parse(f, 'fasta'):
        # do stuff, depending on alphabet

我现在的问题是我不知道我会在.fasta 文件中找到什么样的序列。它可以是蛋白质、DNA 或 RNA，但我必须知道字母表中的字母数量。

有没有办法用 Biopython 从序列中“估计”字母表？我知道一个蛋白质可能只包含字母 ACGT，这就是我想估计字母表的原因。

【问题讨论】：

标签： python bioinformatics biopython fasta

【解决方案1】：

这对于小序列很难做到。例如，序列ACGCGACAGA 可以是 DNA、RNA 或蛋白质序列，因为字母 A、C 和 G 是所有三个字母表共有的。没有其他知识，无法估计哪个是最佳匹配。

以下代码将打印出给定 FASTA 文件中第一条记录所属的所有字母：

from Bio import SeqIO
from Bio.Alphabet.IUPAC import *

alphabets = [extended_protein,  ambiguous_dna, unambiguous_dna, extended_dna, ambiguous_rna, unambiguous_rna]

def validate(seq, alphabet):
    "Checks that a sequence only contains values from an alphabet"
    # inspired by https://www.biostars.org/p/102/
    leftover = set(str(seq).upper()) - set(alphabet.letters)
    return not leftover

with open("example.fasta") as handle:
    first_record = list(SeqIO.parse(handle, "fasta"))[0]
    valid_alphabets = [str(alphabet) for alphabet in alphabets if validate(first_record.seq, alphabet)]
    print("Valid alpahabet(s) for fasta file: {}".format(', '.join(valid_alphabets)))

因此，对于序列ACGCGACAGA，这将打印：

Valid alpahabet(s) for fasta file: ExtendedIUPACProtein(), IUPACAmbiguousDNA(), IUPACUnambiguousDNA(), ExtendedIUPACDNA(), IUPACAmbiguousRNA(), IUPACUnambiguousRNA()

但是对于序列MKQHKAMIVALIVICITAVVAALVTRKDLCEVHIRTGQTEVAVFX，它将打印：

Valid alpahabet(s) for fasta file: ExtendedIUPACProtein()

【讨论】：

这已经解决了我的问题。谢谢你。我知道我不能期待更多，但在大多数情况下，当我有很多序列时，你的代码几乎可以确定。