将核苷酸位置与fasta-fi序列相匹配

Question

我有职位清单：

chr1 1000
chr2 2000
chr3 4000

并希望能够转换这些位置在他们的核苷酸序列给出一个自定义的fasta文件。例如：

^{pr2}$

有没有已经用python编写的工具可以完成这项工作？在

Answer 1

给定FASTA文件chromosomes.fasta：

>chr1
GATTACA
>chr2
ATTACGA
>chr3
GCCAACG

以及位置文件positions.txt：

^{pr2}$

您可以使用以下代码：

from Bio import SeqIO
record_dict = SeqIO.to_dict(SeqIO.parse('chromosomes.fasta', "fasta"))

chromosome_positions = {}
with open('positions.txt') as f:
    for line in f.read().splitlines():
        if line:
            chromosome, position = line.split()
            chromosome_positions[chromosome] = int(position)


for chromosome in chromosome_positions:
    seq = record_dict[chromosome]
    position = chromosome_positions[chromosome]
    base = seq[position]
    print chromosome, position, base

将输出：

chr3 5 C
chr2 4 C
chr1 3 T

注意Python使用zero-based indexing，因此positions.txt中的5位置将给您相应序列中的第六个基。在