通过坐标快速获取人类基因组序列的方法

Question

我想随机获取很多人类基因组片段（超过5亿个）。

这是整个过程的一部分。我有一个来自bowtie的.sam结果文件，里面有1000万个与人类基因组的比对结果。我想把每个查询的读取与.sam文件中它所比对的“参考序列”进行比较。我使用的参考序列是来自UCSC的hg19.fa。所以我需要能够通过.sam文件中的位置信息，从hg19.fa（或染色体文件）中获取相应的序列。

比如，给定位置：chr4:35654-35695，我可以获取到42个碱基的序列：

gtcttccagggtttttatatttttgggttttacacttaagt

到目前为止，我有两个解决方案：

1. 使用python脚本从UCSC DAS服务器获取序列： http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr4:35654,35695
2. 使用python脚本调用“samtools faidx”命令并返回命令输出， 来自帖子： http://seqanswers.com/forums/showthread.php?t=23606&highlight=fetch+genome+coordinate

但是，它们都比较慢。samtools faidx比从DAS服务器获取的速度稍快，但仍然很慢。

那么，有没有什么快速的方法来做到这一点？我有单独的染色体fasta文件和hg19.fa文件。

Answer 1

你可以试试这个叫做python twobitreader的模块：

在命令行输入：python -m twobitreader hg19.2bit < temp.bed

更多信息可以查看这个链接：http://pythonhosted.org/twobitreader/

Answer 2

使用 ucsc 的 twoBitToFa 工具，可以在这个链接找到：http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/

更多信息请查看：http://genome.ucsc.edu/goldenPath/help/twoBit.html