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

Question

如何使用Python高效获取基因组序列？例如，从.fa文件或其他容易获得的格式？我基本上想要一个接口fetch_seq(铬，链，开始，结束)，它将返回序列开始，结束在指定的链上的给定染色体上。

类似地，有没有可编程的python接口来获取phastCons分数？

谢谢。

Answer 1

请看我在Biostar网站上对你问题的回答：

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

对Fasta文件使用SeqIO，您将获得文件中每个项目的记录对象。然后您可以执行以下操作：

region = rec.seq[start:end]

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

Answer 2

从大型人类染色体文件中检索序列数据在内存方面可能效率低下，因此，如果您正在寻找计算效率，可以将序列数据格式化为压缩的二进制字符串，并基于字节位置进行查找。我用perl (可用的here )写了一些例程来做这件事，而python也有同样的pack and unpack routines -所以它是可以做到的，但只有当你在有限的机器上遇到大文件的麻烦时，它才是值得的。否则使用biopython SeqIO

Answer 3

看看biopython，它支持几种基因序列格式。具体来说，它支持FASTA and GenBank files，举几个例子。