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

Question

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

这是整个过程的一部分工作。我有来自蝴蝶结的.sam结果文件，有1000万个人类基因组读取比对。我想将每个查询读取与sam文件中的“它对齐的引用序列”进行比较。我使用的参考序列是UCSC的hg19.fa。所以我需要能够从hg19.fa（或染色体文件）中获取序列，方法是使用sam文件中的位置。在

例如，使用give:chr4:35654-35695，我可以得到42bp序列：

gtcttccagggtttttattttgggttttaacttaagt

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

1. 使用python脚本调用“samtools faidx”命令并返回commnad输出， 发帖人： http://seqanswers.com/forums/showthread.php?t=23606&highlight=fetch+genome+coordinate

但是，他们很慢。samtools faidx比从DAS服务器获取它要快一些，但仍然很慢。在

那么，有没有什么方法可以做到这一点？我有seprate染色体fasta文件和hg19.fa文件。在

Answer 1

您可以尝试使用python twobitreader模块：

python-m双位阅读器hg19.2bit<；临时床在

http://pythonhosted.org/twobitreader/

Answer 2

在http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/中使用ucsctwoBitToFa

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