python - 如何使用Python根据坐标信息检索FASTA序列？

我想根据床文件B.bed获取序列，该文件包含序列的坐标信息，方法是将坐标与A.fasta的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')

快餐

>chr16:13361561-13361573
TAGTGGGTCAGAC
>chr6_apd_hap1:2165669-2165681
AGATGAGTCATCA
>chr10:112612173-112612185
AAGTGTGTCAGCT

床

chr6_apd_hap1   2165668 2165681
chr10   112612172   112612185

C.fasta的预期输出

>chr6_apd_hap1:2165669-2165681
AGATGAGTCATCA
>chr10:112612173-112612185
AAGTGTGTCAGCT

但是我收到以下错误：

long_seq_record = records[name]
KeyError: 'chr6_apd_hap1'

最佳答案

在您的代码中，positions是一个defaultdict，它具有BED文件中的名称作为键：

>>> 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='')

关于python - 如何使用Python根据坐标信息检索FASTA序列？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/31285738/