在Python中有效获取基因组序列？

Question

如何使用 Python 高效获取基因组序列？例如，从 .fa 文件或其他一些容易获得的格式？我基本上想要一个接口 fetch_seq(chrom,strand,start,end) 它将返回指定链上给定染色体上的序列 [start,end]。

类似地，是否有一个用于获取 phastCons 分数的编程 Python 接口？

谢谢。

Answer 1

请参阅我在 Biostar 对您问题的回答：

http://biostar.stackexchange.com/questions/1639/getting-genomic-sequences-and-phastcons-scores-using-python-from-ensembl-ucsc

将 SeqIO 与 Fasta 文件结合使用，您将获得文件中每个项目的记录对象。然后你可以这样做：

region = rec.seq[start:end]

拉出切片。使用标准库的好处是您不必担心原始 fasta 文件中的换行符。

Answer 2

从大型人类染色体文件中检索序列数据在内存方面可能效率低下，因此如果您正在寻求计算效率，您可以将序列数据格式化为打包二进制字符串并根据字节位置进行查找。我用 perl 编写了例程来执行此操作（可用 这里 ），和python有相同的 打包和解包例程 - 所以这是可以做到的，但只有当您在有限的机器上遇到大文件的麻烦时才值得。否则使用biopython SeqIO

Answer 3

看一眼 生物蟒蛇, ，它支持多种基因序列格式。具体来说，它支持 FASTA 和 GenBank 文件, ，举几个例子。

Answer 4

皮法斯塔 是您正在寻找的模块。从描述来看

快速、内存高效、Pythonic（和命令行）访问 fasta 序列文件

https://github.com/brentp/pyfasta