Python:使用Bed文件从FASTA文件中提取DNA序列

Question

我可以知道如何从fasta文件中提取dna序列吗？我试过床上工具和工具。贝德工具公司做得不错，但我的一些文件返回了“警告:在fasta文件中没有发现染色体”，但事实是床文件中的染色体名与fasta完全相同。我正在寻找其他可供选择的选项，python可以为我完成此任务。

床铺档案：

chr1 1:117223140-117223856 3 7

chr1 1:117223140-117223856 5 9

Fasta档案：

chr1 1:117223140-117223856

CGCGTGGGCTAGGGGCTAGCCCC

期望产出：

chr1 1:117223140-117223856

CGTGG

chr1 1:117223140-117223856

TGGGC

Answer 1

BioPython是您想要使用的：

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('positions.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('sequences.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(short_seq, alphabet), id=name, description='')
        short_seq_records.append(short_seq_record)

# write to file
with open('output.fasta', 'w') as f:
    SeqIO.write(short_seq_records, f, 'fasta')

Answer 2

您的床文件需要标签分隔的床工具使用它。用选项卡替换冒号、破折号和空格。

BedTools文档页面说：“床工具要求所有床输入文件(以及从stdin接收的输入)都是制表符分隔的。”BedTools。