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

Question

我该如何用Python高效地获取基因组序列呢？比如说，从一个.fa文件或者其他容易获取的格式中。我基本上想要一个接口，叫做fetch_seq(chrom, strand, start, end)，它能返回在指定染色体和方向上，从开始位置到结束位置的序列。

类似地，有没有什么Python接口可以用来获取phastCons分数呢？

谢谢。

Answer 1

看看这个 biopython，它可以支持几种基因序列的格式。特别是，它支持 FASTA 和 GenBank 文件，这只是其中的两个例子。

Answer 2

从大型人类染色体文件中获取序列数据可能会占用很多内存，所以如果你想提高计算效率，可以把序列数据格式化成一种紧凑的二进制字符串，然后根据字节位置来查找。我写了一些程序来实现这个功能，使用的是perl语言（可以在这里找到），而python也有类似的打包和解包功能 - 所以这是可以做到的，但只有在你遇到大型文件在有限的机器上处理困难时才值得这样做。否则可以使用biopython的SeqIO模块。

Answer 3

你可以看看我在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文件中的换行符。