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

Question

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

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

谢谢。

Answer 1

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

Answer 2

在 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 3

看看biopython，它支持多种基因序列格式。具体来说，它支持FASTA 和 GenBank 文件，仅举几例。

Answer 4

pyfasta是您正在寻找的模块。从描述

快速、节省内存、pythonic（和命令行）访问 fasta 序列文件

https://github.com/brentp/pyfasta