如何使用 biopython 查询 fasta 文件中的区域？答案

【问题标题】：How to query a region in a fasta file using biopython?如何使用 biopython 查询 fasta 文件中的区域？
【发布时间】：2021-01-22 08:38:33
【问题描述】：

我有一个包含一些参考基因组的 fasta 文件。在给定染色体、开始和结束索引的情况下，我想将参考核苷酸作为字符串获取。

我正在寻找一个在代码中看起来像这样的函数：

from Bio import SeqIO
p = '/path/to/refernce.fa'
seqs = SeqIO.parse(p.open(), 'fasta')
string = seqs.query(id='chr7', start=10042, end=10252)

字符串应该是这样的：'GGCTACGAACT...'

我发现的只是如何迭代 seq，以及如何从 NCBI 中提取数据，这不是我想要的。在 biopython 中执行此操作的正确方法是什么？

【问题讨论】：

这是否适用于 biopython next(r for r in seqs if r.id == "chr7").seq[10042:10252]？注意它可能是 end+1
可能会，但需要一段时间，下面的解决方案非常好:) 非常感谢您的回复想法！

【解决方案1】：

AFAIK，biopython 目前没有此功能。对于使用索引的随机查找（请参阅samtools faidx），您可能需要pysam 或pyfaidx。这是一个使用 pysam.FastaFile 类的示例，它允许您快速“获取”区域中的序列：

import pysam
ref = pysam.FastaFile('/path/to/reference.fa')
seq = ref.fetch('chr7', 10042, 10252)
print(seq)

或者使用pyfaidx和'get_seq'方法：

from pyfaidx import Fasta
ref = Fasta('/path/to/reference.fa')
seq = ref.get_seq('chr7', 10042, 10252)
print(seq)

【讨论】：