我正在寻找一种.fasta
在Biopython中读取文件的方法,如果我们处理的是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,这就是为什么我想估计字母表.
对于小序列来说这很难做到.例如,序列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()