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

Question

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

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

例如通过给出:chr4:35654-35695，我可以获得 42bp 序列:

gtcttccagggttttttatattttggggtttacacttaagt

到目前为止，我有 2 个解决方案: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 服务器获取它快一点，但仍然很慢。

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

Answer 1

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

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