一种利用坐标快速获取人类基因组序列的方法

Question

我想随机获得很多人类基因组片段(其中5亿多个)。

这是整个过程的部分工作。我有.sam的结果文件从保龄球，与1000万人类基因组读对齐。我希望将每个查询读取与sam文件中的“引用序列”进行比较。我使用的参考序列是来自UCSC的hg19.fa。因此，我需要能够从hg19.fa (或染色体文件)中获取序列，使用sam文件中的位置。

例如，通过给予:chr4 4:35654-35695，我可以得到42 get的序列：

gtcttccagggtttttatatttttgggttttacacttaagt

到目前为止，我有两个解决方案: 1.从UCSC服务器获取序列的python脚本：http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr4:35654,35695

1. 使用python脚本调用‘samtools’命令并返回commnad输出，从post：http://seqanswers.com/forums/showthread.php?t=23606&highlight=fetch+genome+coordinate

但是，他们很慢。samtools比从DAS服务器获得它要快一些，但速度仍然很慢。

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

Answer 1

在ucsc twoBitToFa中使用64/

另见http://genome.ucsc.edu/goldenPath/help/twoBit.html

Answer 2

您可以尝试模块：

python -m双读器hg19.2位< temp.bed

http://pythonhosted.org/twobitreader/