a fast way to get human genome sequence by coordinate

[xiangwulu]

Hi all,

I want to get a lot human genome fragments (more than 500 million of them) randomly.


This is a partial work of the whole process. I have .sam result file from bowtie, with 10 million human genome reads alignment. I want to compare each query reads with the 'reference sequence it aligned to' from the sam file. The reference sequence I used is hg19.fa from UCSC. So I need to be able to get the sequence from hg19.fa (or chromosome files) by using the location in the sam file.


e.g. with giving: chr4:35654-35695, i could get 42bp sequences:
gtcttccagggtttttatatttttgggttttacacttaagt

so far, i had 2 solutions:
1. python script to fetch sequences from UCSC DAS server:

http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr4:35654,35695

2. using python script call ''samtools faidx'' command and return commnad output,
from post:
http://seqanswers.com/forums/showthr...ome+coordinate

but, they are slow. samtools faidx is bit faster than getting it from DAS server, but still slow.

so, is there any FAST way to do this? i have the seprate chromosome fasta files, and hg19.fa file.

[swbarnes2]

Doesn't the MD line tell you how the read differs from the genome? Wouldn't it be easier to use that and the sequence from the SAM line, so you don't have to look at another file at all?

[westerman]

How fast do you want it to be? 'samtools faidx', in my experience, is rather speedy but you are asking for a lot of data. As a test, using bash I tried fetching 100,000 genome coordinates (randomly generated asking for lengths from 50 to 100 bases on all chromosomes) via 'samtools faidx'. It took about 7 minutes wall-clock time and generated a 10 MB file. At that rate your desired 500 million would take ... hum, around 25 days and generate a 50GB file.

If I do the 100,000 coordinates 100 at a time (i.e., 100 regions on the 'samtools faidx' command line thus avoiding having to start up samtools 100,000 times as opposed to only 1000 times) then the time remains the same. So I suspect that the delay is either samtool's efficiency or simply having to read the largish human genome.

You might be able to cut that time down if you first separated your input into chromosomes and then did per-chromosome 'samtools faidx' in parallel -- assuming that your disks are fast enough. And/or you can use a RAM-disk to hold the genome file and the *.fai file.

I doubt if doing a web fetch on 500 million fragments will be fast nor will make you popular with UCSC.

jm's solution could also be feasible. No matter how you slice it taking 500 million of anything takes time.

[westerman]

swbarne's comment (use the MD line) was also something I was thinking about. But it takes more programming than using samtools directly.

[vivek_]

System calls from within a scripting language will always add a lot of latency. Since you seem interested in using Python, I came across this module yesterday which might be applicable here

https://github.com/mdshw5/pyfaidx

Samtools provides a function "faidx" (FAsta InDeX), which creates a small flat index file ".fai" allowing for fast random access to any subsequence in the indexed fasta, while loading a minimal amount of the file in to memory.

Pyfaidx provides an interface for creating and using this index for fast random access of DNA subsequences from huge fasta files in a "pythonic" manner. Indexing speed is comparable to samtools, and in some cases sequence retrieval is much faster (benchmark).

[dpryan]

If you need it to be quick then just read the genome into memory once and then iterate over things. That's how those of us who have written methylation callers (where you have to do this exact process) do things, since it gives the best performance.

Anyway, as the others said, it's easier to just parse the MD string if it exists and is valid.

[swbarnes2]

swbarne's comment (use the MD line) was also something I was thinking about. But it takes more programming than using samtools directly.

In the 25 days it would take to use samtools faidx hundreds of millions of times, you could learn how to program well enough to parse the SAM entry.

[xiangwulu]

Doesn't the MD line tell you how the read differs from the genome? Wouldn't it be easier to use that and the sequence from the SAM line, so you don't have to look at another file at all?

The MD tag in sam files tell something, but not enough. Also, some tools does not output MD tag in sam file, like blat, bfast.

[mdshw5]

System calls from within a scripting language will always add a lot of latency. Since you seem interested in using Python, I came across this module yesterday which might be applicable here

https://github.com/mdshw5/pyfaidx

This is exactly right. System calls from a scripting language will be slow, and pyfaidx is quite a bit better than calls to samtools or pysam (see Table 1 in my pre-print for pyfaidx). However, dpryan's suggestion of using raw strings (or something close, like Biopython's Seq.IO) will always be the fastest option, at the expense of keeping track of your own substring indexes. If "fast enough" is okay, I would suggest using pyfaidx as it's extensively tested for correctness and has a reasonable API for working with your existing indexed fasta files efficiently.

However, I may be biased.