Also, it simplifies random access of the base stream directly from the file. You don't have to keep track of where you really are based on all the previous insertions.

The idea that .2bit files are more compressed is probably in relationship to the earlier .nib format which uses 4 bits per base. Clearly the frequency of N and M bases is such that the 2bit format is much more compact.

I think that a virtual api to a 2bit file would be useful. So you could do something like:

my $twobit = new TwoBitFile('some/path.2bit');
# returns 1000 bases starting from base 40,000,000:
my $bases = $twobit->range(40_000_000, 1000);
[download]

Update: Being able to treat the 2bit file as a virtual string with the ability to run regular expressions on it would also be really cool. Is that something we can do in Perl6?

pc88mxer,
With regards to the N blocks, your hypothesis sounds valid. With regards to the virtual API, I was actually thinking of turning it into an iterator. Calling the iterator without any parameters would fetch the next X bytes where X is some default. You could then allow for named parameters to fetch forward by Y, rewind, select a range, etc. It should have gone into the TODO but wasn't really necessary.

With regards to Perl 6 and applying a regular expressions on streams - the answer is yes unless things have changed since last I looked.

Cheers - L~R

my @CONV = glob( "{T,C,A,G}" x 4 );
my $dna = join "", @CONV[ unpack "C*", $raw ];
[download]

## takes 16 bits (= 8 bases = unsigned short) at a time
my @CONV = glob( "{T,C,A,G}" x 8 );
my $dna = join "", @CONV[ unpack "S*", $raw ];
[download]

Another cute trick I can think of is that you can do some bit-twiddling to implement the M-blocks (apparently lowercasing a range of characters). In ASCII, you can toggle the case of an alphabetic character by bitwise-XOR'ing it with the space character. So I think you can rewrite:

substr($dna, $_, $mblock{$_}, lc(substr($dna, $_, $mblock{$_})))
[download]

substr($dna, $_, $mblock{$_}) ^= (" " x $mblock{$_});
[download]

Of course you could always fix M,N blocks on-the-fly, as you are unpacking them from $raw, but that would require some more work. Since I'm typing one-handed these days and it takes me forever, I think I will pass on playing with some code that does that! ;)

Your idea of looking up the meanings of the sequences byte per byte is brilliant. I do have some doubts about using glob for it... but it even appears to do the right thing on ActivePerl on Windows. Still, I'm wondering if this is not just pure luck.

A reliable way to do it would be to generate a list of integers, in this case from 0 to 255, and convert each to a string using base 4 — admittedly, I don't know how to best do it in Perl... As a second step, I'd convert the digits '0' .. '3' to the letters, for example with

tr/0123/TCAG/
[download]

Anyway, you say

at the expense of having a lookup table for all 256 values

WTF? What expense is that? A few k of memory? Seriously, if the proper way to generate the array of meanings is too expensive, I'd just generate it once at startup, and store it in memory.

You could also experiment with different trade-offs on lookup table sizes:

Yes, but in that case, the lookup table gets much bigger: 64k entries of 8 letters each, that is 256k of text plus overhead of the array. Ouch. I don't think it will be much faster, so I don't think it's worth it.

For some reason, I had byte-order issues doing this.

Of course you have. You used a machine dependent byte ordering. You should either use 'n' or 'v' as the basic unpack template (probably 'n', for Big Endian), which luckily appears to produce unsigned integers, too.

I have some doubts about using unpack "C*", $raw to convert the byte sequence into numbers. Ouch. That sequence can be millions or even billions of bytes long, and that is a very long list. I think it's better to convert the $raw string either in short sequences of, say, a few k each (the compromise is in loop count vs. memory usage per loop),

my $dna = '';
use constant CHUNKSIZE => 2048;
for (my $offset = 0; $offset < length($raw); $offset += CHUNKSIZE) {
   $dna.= join '', @CONV[ unpack 'C*', substr $raw, $offset, CHUNKSIZE
+ ];
}
[download]

or maybe even byte per byte with s///:

s/(.)/$CONV[ord $1]/sge
[download]

but I doubt this will be the fastest way. It'll be as memory cheap as possible, that is true.

Finally: don't forget to cut off the junk at the end of the sequence, making the length the same as the number of entries there were expected according to the record header.

The first thing I noticed is that you are detecting the byte order for the header, but then ignoring it and using the platform specific 'l' template then on. That's wrong in two ways:

1. Using 'l' template, which is for signed 32-bit numbers could cause problems.

   You should probably at least use 'L'. But I think you ought to be using 'N' or 'V'.

2. My interpretation on the spec you linked is that the entire file may need byte order reversal?

   If you've already successfully parsed a 2bit file, then you may have just lucked out on using the same platform/byteordering as it was created on.

   My instinct would be to check the sig and save a template char of 'N' or 'V' as appropriate and use that for all future unpacking. Has a knock-on advantage also. See later.

The slowest bit of the process seems likely to be

    $dna = join '', 
           map $conv{$_},
           unpack('(A2)*', unpack("B" . $bytes * BITS_PER_BYTE , $raw)
+);
[download]

You should be able to save a bit of time by building a larger lookup table:

my %bitMap  =  ('00' => 'T', '01' => 'C', '10' => 'A', '11' => 'G');
my @byteMap = map{ 
    join '', map $bitMap{ $_ }, unpack '(A2)4', unpack 'B8', chr 
} 0 .. 255;
[download]

Now you can convert each byte (4 bases) in the packed DNA to it's ascii with a single array lookup rather than 4 hash loopkups:

## Omit the braces add a comma for a negligable further improvement
my $DNA = join '', map{ $byteMap[ $_ ] } unpack 'C*', $raw;
[download]

Also, build your lookups at compile time not over and over at runtime as now.

Of course, you can take that idea a little further and do two bytes at a time:

my %bitMap  =  ('00' => 'T', '01' => 'C', '10' => 'A', '11' => 'G');
my @byteMap = map{ 
    join '', map $bitMap{ $_ }, unpack '(A2)4', unpack 'B8', chr 
} 0 .. 255;
my @wordMap = map {
    $byteMap[ $_ >> 8 ] . $byteMap[ $_ | 0xff ]
} 0 .. 65535;

...

my $DNA = join '', map{ $wordMap[ $_ ] } unpack 'n*', $raw;
[download]

Which ought to be close to an order of magnitude faster than your current method. 1 array lookup -v- 8 hash lookups; 8 times less lower loop overhead.

Of course, the 'v' template will be byte-order specific. But, if when you determine the byte order from reading the signature, you store a template of 'N' or 'V' for your 32-bit field processing, then you can just lc that template to obtain the unsigned short template.

Also, how sure are you of your conversion table? Are you certain that you should be using 'B' and not 'b'.

I might have tested some of this, but it would take me the best part of two days to download the "sample" .2bit file you linked and a quick google didn't locate any others.

Final thought. If you built an ordered array of offsets, as well as the named index, when processing the toc, you could provide access by position as well as access by name. (Your ordered hash module might work for this also :)

ftp://hgdownload.cse.ucsc.edu/gbdb/ce4/ce4.2bit is much smaller, 24.6MB

Nice example++. Thanks for posting the code.

The 2bit format is used by the BLAT program for aligning DNA sequences. The BLAT suite of programs contains two utilities that might be of interest:

• faToTwoBit – convert Fasta format sequence files to a dense randomly-accessable .2bit format
• twoBitToFa – convert from the .2bit format back to fasta

The twoBitToFa documentation is as follows:

twoBitToFa - Convert all or part of .2bit file to fasta

usage:
   twoBitToFa input.2bit output.fa

options:
   -seq=name - restrict this to just one sequence
   -start=X  - start at given position in sequence (zero-based)
   -end=X - end at given position in sequence (non-inclusive)
[download]

Once the sequence is converted back to fasta format it can be easily manipulated using a variety of tools (Perl, BioPerl, sequence utility programs, etc).

So yes, a utility does exist that meets at least some of your listed requirements. OTOH, another tool in the toolbox, especially if it is pure Perl and/or a snazzy wrapper around a command-line utility, is always welcome. :-)

    # Parse header
    my ($sig, $ver, $cnt, $reserved) = unpack('l4', $raw);

    # TODO: validate (signature, reverse byte order, version)
    return {SIG => $sig, VER => $ver, CNT => $cnt, RSV => $reserved};
[download]

All fields are 32 bits unless noted. If the signature value is not as given, the reader program should byte-swap the signature and check if the swapped version matches. If so, all multiple-byte entities in the file will have to be byte-swapped. This enables these binary files to be used unchanged on different architectures.

for my $template ("N", "V") {
    # Parse header
    my ($sig, $ver, $cnt, $reserved) = unpack($template.'4', $raw);
    if($sig==0x1A412743) {
        return {unpack => $template, VER => $ver, CNT => $cnt};
    }
}
# no match: not a .2bit header
return undef;
[download]

Now, that wasn't so hard, was it?

For the rest... See blokhead's node — and my reply with my remarks. I agree with his approach as it's probably as fast as you can get in pure Perl. I'd stick to decoding whole bytes in one go. Be careful about memory problems, though: these strings can be very long — and your own approach is even worse, as it uses strings double that size, for decoding, and you're copying it around some more, taking up even more temporary space. Ouch.

And using more memory usually means (much) slower — if you're not simply running out of memory, or you have to close some other programs to be able to run yours. None of them good things.

I'm sorry to tell you that if you compare the output generated from your code and the one from the original twobittofa program, they don't match. So something's wrong...

Anonymous Monk,
No, nothing is wrong but thank you for your response. The twobittofa program, which I was not attempting to duplicate converts a .2bit file into a fasta file as I understand it. The intent of my project was not to duplicate that program - if that's what you need, then that's what you should use.

Cheers - L~R