Python：使用 Bed 文件从 FASTA 文件中提取 DNA 序列答案

【问题标题】：Python: Extract DNA sequence from FASTA file using Bed filePython：使用 Bed 文件从 FASTA 文件中提取 DNA 序列
【发布时间】：2015-05-28 10:12:22
【问题描述】：

我可以知道如何从 fasta 文件中提取 dna 序列吗？我尝试了床具和 samtools。 Bedtools getfasta 做得很好，但对于我的一些文件返回“警告：在 fasta 文件中找不到染色体”，但事实是 bed 文件中的染色体名称和 fasta 完全相同。我正在寻找 python 可以为我完成此任务的其他替代方法。

床文件：
chr1:117223140-117223856 3 7
chr1:117223140-117223856 5 9

快速文件：
>chr1:117223140-117223856
CGCGTGGGCTAGGGGCTAGCCCC

所需输出：
>chr1:117223140-117223856
CGTGG
>chr1:117223140-117223856
TGGGC

【问题讨论】：

您想要的输出似乎有误。染色体从 0 开始编号——从 BED format description 开始——因此从 cr1:117223140-117223856 中提取编号为 3 到 7 的染色体将得到 GTGGG，从中提取编号为 5 到 9 的将是 GGGCT。

标签： python bioinformatics biopython fasta bed

【解决方案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')

【讨论】：

谢谢。这很有帮助。如何确保输出顺序遵循bed文件？

【解决方案2】：

尝试，用：

from Bio import SeqIO

#I use RAM, and to store fasta in dictionary
parser = SeqIO.parse(open("input.fasta")
dict_fasta = dict([(seq.id, seq) for seq in parser, "fasta")])

output = open("output.fasta", "w")
for line in open("input.bed"):
  id, begin, end = line.split()
  if id in dict_fasta:
    #[int(begin)-1:int(end)] if the first base in a chromosome is numbered 1
    #[int(begin):int(end)+1] if the first base in a chromosome is numbered 0
    output.write(dict_fasta[id][int(begin)-1:int(end)].format("fasta"))
  else:
    print id + " don't found"

output.close()

你知道，染色体的第一个碱基编号为 1：

>chr1:117223140-117223856 CGTGG >chr1:117223140-117223856 TGGC

你知道，染色体的第一个碱基编号为 0：

>chr1:117223140-117223856 GTGGG >chr1:117223140-117223856 GGGCT

【讨论】：

【解决方案3】：

您的床文件需要用制表符分隔，床具才能使用它。用制表符替换冒号、破折号和空格。

BedTools 文档页面显示“bedtools 要求所有 BED 输入文件（以及从标准输入接收的输入）都是制表符分隔的。” BedTools.

【讨论】：