将核苷酸位置与 fasta 文件中的序列匹配

Question

我有职位列表：

chr1 1000
chr2 2000
chr3 4000

并希望能够转换其核苷酸序列中的那些位置，从而提供自定义的 fasta 文件。如：

chr1 1000 A
chr2 2000 T
chr3 4000 G

python 中是否有任何已经编写好的工具可以完成这项工作？

Answer 1

给定 FASTA 文件 chromosomes.fasta:

>chr1
GATTACA
>chr2
ATTACGA
>chr3
GCCAACG

和位置文件 positions.txt:

chr1 3

chr2 4

chr3 5

您可以使用以下代码：

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会给你相应序列中的第6个碱基。

将核苷酸位置与 fasta 文件中的序列匹配

Match nucleotide position to sequence from fasta file

python

bioinformatics

fasta

biopython