如何使用Python根据坐标信息检索FASTA序列？
如何使用Python根据坐标信息检索FASTA序列？
提问于 2015-07-08 07:08:50
我希望根据包含序列坐标信息的床文件B.bed，通过将坐标匹配到fasta文件(即A.fasta )获得序列序列，并根据B.bed文件检索相应的序列。Fasta文件是具有序列ID的正常文件，其次是其核苷酸ATCG。但我得到了下面的关键错误。有谁可以帮我？
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from collections import defaultdict
# read names and postions from bed file
positions = defaultdict(list)
with open('B.bed') as f:
for line in f:
name, start, stop = line.split()
positions[name].append((int(start), int(stop)))
# parse faste file and turn into dictionary
records = SeqIO.to_dict(SeqIO.parse(open('A.fasta'), 'fasta'))
# search for short sequences
short_seq_records = []
for name in positions:
for (start, stop) in positions[name]:
long_seq_record = records[name]
long_seq = long_seq_record.seq
alphabet = long_seq.alphabet
short_seq = str(long_seq)[start-1:stop]
short_seq_record =SeqRecord(Seq(name))+ '\t'+str(start)+'\t'+str(stop)+'\t' + SeqRecord(Seq(short_seq, alphabet))
short_seq_records.append(short_seq_record)
# write to file
with open('C.fasta', 'w') as f:
SeqIO.write(short_seq_records,f, 'fasta')A.fasta
>chr16:13361561-13361573
TAGTGGGTCAGAC
>chr6_apd_hap1:2165669-2165681
AGATGAGTCATCA
>chr10:112612173-112612185
AAGTGTGTCAGCTB.bed
chr6_apd_hap1 2165668 2165681
chr10 112612172 112612185C.fasta的预期输出
>chr6_apd_hap1:2165669-2165681
AGATGAGTCATCA
>chr10:112612173-112612185
AAGTGTGTCAGCT但我收到了以下错误：
long_seq_record = records[name]
KeyError: 'chr6_apd_hap1'回答 1
Stack Overflow用户
回答已采纳
发布于 2015-07-08 08:28:36
在您的代码中，positions是一个defaultdict，它以床文件中的名称作为键：
>>> print positions.keys()
['chr10', 'chr6_apd_hap1']records是一种字典，它以FASTA文件的头作为键，减去开头的>，但它们仍然包括冒号和染色体上的位置：
>>> print records.keys()
['chr16:13361561-13361573', 'chr6_apd_hap1:2165669-2165681', 'chr10:112612173-112612185']因此，您首先需要转换您的records键来释放额外的信息，以便您可以使用positions键来检索它们。您可以在创建records字典之后添加以下行来完成此操作：
records = {key.split(':')[0]: value for (key, value) in records.iteritems()}而且，您当前构建short_seq_record的方式并不真正有效。更换线路
short_seq = str(long_seq)[start-1:stop]
short_seq_record =SeqRecord(Seq(name))+ '\t'+str(start)+'\t'+str(stop)+'\t' + SeqRecord(Seq(short_seq, alphabet))通过以下方式：
short_seq = Seq(str(long_seq)[start-1:stop], alphabet)
short_seq_id = '{0}\t{1}\t{2}'.format(name, start, stop)
short_seq_record = SeqRecord(short_seq, id=short_seq_id, description='')页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/31285738
复制相关文章
相似问题
腾讯云开发者
