我想随机获得很多人类基因组片段（超过 5 亿个）。
这是整个过程的部分工作。我有来自 bowtie 的 .sam 结果文件，有 1000 万个人类基因组读取对齐。我想将每个查询读取与 sam 文件中的“它对齐的参考序列”进行比较。我使用的参考序列是来自 UCSC 的 hg19.fa。所以我需要能够通过使用 sam 文件中的位置从 hg19.fa （或染色体文件）中获取序列。
例如，给予：chr4:35654-35695，我可以获得 42bp 序列：
gtcttccagggtttttattatttttgggttttacacttaagt
到目前为止，我有 2 个解决方案： 1. 从 UCSC DAS 服务器获取序列的 python 脚本：http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr4: 35654,35695
- 使用 python 脚本调用“samtools faidx”命令并返回 commnad 输出，来自帖子： http ://seqanswers.com/forums/showthread.php?t=23606&highlight=fetch+genome+coordinate
但是，它们很慢。samtools faidx 比从 DAS 服务器获取它要快一些，但仍然很慢。
那么，有什么快速的方法可以做到这一点吗？我有单独的染色体 fasta 文件和 hg19.fa 文件。