Dear Assaf and Fed,
Thanks for your suggestion. It has been incorporated in the latest changeset (#2488).
Guru Galaxy team.
On Jul 17, 2009, at 1:08 PM, Assaf Gordon wrote:
Hello,
We've encountered a small issue when using the extract-genomic- sequences tool with a TwoBit file and lot's of sequences: memory and CPU usage increases unnecessarily (and processing time, too).
Background: a TwoBit (2bit) file is a binary representation of a FASTA file. Unlike nbit files, TwoBit files contain multiple sequences per file.
For example, The Drosophila Virilis 'genome' contains 13530 scaffolds. Converting the droVir3 FASTA file into a 2bit file creates a 52MB file.
Details: When trying to extract sequences (using extract_genomic_dna.py), Each time a new chromosome/scaffold is encountered, the *entire* 2bit file is loaded and stored. This is unnecessary - the same 2bit file is used for all the scaffolds. (Imagine a BED file with 1000 different scaffolds - memory usage goes to 1000 * 52MB... ).
The offending code is (extract_genomic_dna.py) :
elif twobit_path and os.path.exists( twobit_path ): if chrom in twobits: t = twobits[chrom] else: twobits[chrom] = t = bx.seq.twobit.TwoBitFile( file( twobit_path ) ) try: sequence = t[chrom][start:end]
the hash 'twobits' contains multiple copies of the same 2bit file, over and over again (one per chromosome/scaffold), but only one chromosome/scaffold is read from it.
An improvement would be to load the file only once, and each time take the chromosome/scaffold from the already-loaded hash.
Here's a possible patch:
--- extract_genomic_dna.py 2009-07-17 12:55:53.000000000 -0400 +++ extract_genomic_dna2.py 2009-07-17 12:55:18.000000000 -0400 @@ -84,6 +84,7 @@ def __main__(): fout = open( output_filename, "w" ) warnings = [] warning = ''
twobitfile = None
for i, line in enumerate( open( input_filename ) ): line = line.rstrip( '\r\n' )
@@ -132,12 +133,10 @@ def __main__(): invalid_line = line continue elif twobit_path and os.path.exists( twobit_path ):
if chrom in twobits:t = twobits[chrom]else:twobits[chrom] = t = bx.seq.twobit.TwoBitFile(file( twobit_path ) )
if not twobitfile:twobitfile = bx.seq.twobit.TwoBitFile( file(twobit_path ) ) try:
sequence = t[chrom][start:end]
sequence = twobitfile[chrom][start:end] except: warning = "Unable to fetch the sequence from'%d' to '%d' for build '%s'. " %( start, end-start, dbkey ) warnings.append( warning ) ========================
Technical note: This code and patch refers to an older revision (before changeset 2457:8dfe971fcc27), but I think the new code still behaves this way.
Kudos to Ted Karginov who had to learn python to find and fix this bug.
-Gordon. _______________________________________________ galaxy-dev mailing list galaxy-dev@bx.psu.edu http://mail.bx.psu.edu/cgi-bin/mailman/listinfo/galaxy-dev
Guruprasad Ananda Graduate Student Bioinformatics and Genomics The Pennsylvania State University