TwoBitFile

TwoBitFile

• TwoBitFile is a python module dealing with huge genome files in 2bit format.
• This class is simplified version of public available TwoBitReader.

Definition of TwoBitFile:

class TwoBitFile(dict):
    ''' This class is used to deal with huge genome files in 2bit format. '''
   7 lines:     def __init__(self,infile):-----------------------------------------------------------------------------------------------
  12 lines:     def _loadheader(self):---------------------------------------------------------------------------------------------------
  20 lines:     def _loadindex(self):----------------------------------------------------------------------------------------------------
  11 lines:     def chromSizes(self):----------------------------------------------------------------------------------------------------
   3 lines:     def __del__(self):-------------------------------------------------------------------------------------------------------

class TwoBitSeq(object):
    ''' TwoBit sequence. '''
  37 lines:     def __init__(self, fh, offset, fsize, byteswapped = False):--------------------------------------------------------------
   3 lines:     def __len__(self):-------------------------------------------------------------------------------------------------------
   3 lines:     def __str__(self):-------------------------------------------------------------------------------------------------------
  60 lines:     def __getslice__(self, _min, _max = None):-------------------------------------------------------------------------------

1. Example get sequence sizes from 2bit file

> from ngslib import TwoBitFile
> g = TwoBitFile("test/test.2bit")
> chroms = g.chromSizes()
> chroms

Output:

{'chr4': 250L, 'chr3': 320L, 'chr2': 240L, 'chr1': 240L}

2. Example get slice of sequences

> chr1 = g['chr1']
> chr1[100:200]

Output:

'AATATGAAGTTCTTTAGCATAACAAGGATCTGCCTTTGTAAAAGAAaaagaaagaaagagcgaaagaaagaaaAGAACTGAGGACAGCATTCTTTTCTCT'

Extract the whole chromosome. Note: this will take a long time if the chromosome is huge.

> str(chr1) # the same as chr1[:]

Output:

'CTCTTACGTTTTATTCCCTCTTTATCTCAGCTTAGATCAGGGTAAACTTTCAGAAAGCCTTACTGTTGCATTTTGTTAGTTTCTGTTTTCCTCAACAACT
AATATGAAGTTCTTTAGCATAACAAGGATCTGCCTTTGTAAAAGAAaaagaaagaaagagcgaaagaaagaaaAGAACTGAGGACAGCATTCTTTTCTCTC
ACCATTGTCAGTGGGTAGGCAAATGCTGTGTATACCTAA'

Data structure of 2bit file

Step 1: Read header of 2bit file.

fh = open("test/test.2bit",'rb')
header = array(LONG)
header.fromfile(fh, 4)
(signature, version, sequence_count, reserved) = header
if signature == 0x1A412743:
    print "version is {0}, sequence count is {1} and reserved value is {2}".format(version, sequence_count, reserved)

fh.seek(16)
offset1=0 # record the offset for the last chromosome
for i in range(sequence_count):
    name_size = array('B')
    name_size.fromfile(fh,1)
    print "name size is",name_size[0],
    name = array('c')
    name.fromfile(fh,name_size[0])
    print ", name is", name.tostring(),
    offset = array(LONG)
    offset.fromfile(fh, 1)
    offset1 = offset[0]
    print ", the offset is", offset1

Output:

version is 0, sequence count is 4 and reserved value is 0
name size is 4 , name is chr1 , the offset is 52
name size is 4 , name is chr2 , the offset is 144
name size is 4 , name is chr3 , the offset is 252
name size is 4 , name is chr4 , the offset is 356

Step 2: Read header of each chromosome.

# go to the last chrom directly
fh.seek(offset1) 

# read header
header = array(LONG)
header.fromfile(fh,2)
dnasize,blockcount = header
print "DNA size is", dnasize
print "Block count is", blockcount

# read blocks
bstarts=array(LONG)
bsizes=array(LONG)
bstarts.fromfile(fh, blockcount)
bsizes.fromfile(fh, blockcount)
for i in range(blockcount):
    print bstarts[i],bsizes[i]

# read masks
maskcount = array(LONG)
maskcount.fromfile(fh, 1)
maskcount = maskcount[0]
print "Mask count is", maskcount
mstarts=array(LONG)
msizes=array(LONG)
mstarts.fromfile(fh, maskcount)
msizes.fromfile(fh, maskcount)
for i in range(maskcount):
    print mstarts[i],msizes[i]

# record the start offset of sequence
fh.read(4)
offset2 = fh.tell()

Output:

DNA size is 250
Block count is 2
17 6
75 5
Mask count is 1
172 78

Step 3: Extract sequence by slice

###########################################
# 1 byte = 4 bases
# _min = 10 is in the 3rd block
# _max = 100 is in the 25th block
# in total we have 23 blocks to read
###########################################
_min=10
_max=100
start_block = _min / 4
end_block = (_max + 3) / 4
blocks_to_read = end_block - start_block
print "start block:", start_block, ", end block:", end_block
print blocks_to_read, "blocks to read."

# read blocks
fh.seek(offset2 + start_block)
fourbyte_dna = array('H')
fourbyte_dna.fromfile(fh,blocks_to_read/2)
last_byte = array('B')
if blocks_to_read % 2 == 1:
    last_byte.fromfile(fh,1)
seq = array('c', "".join([twoBytesTable[i] for i in fourbyte_dna]) + (len(last_byte) and byteTable[last_byte[0]] or ""))
print "Sequence in bytes:\n",seq.tostring()

# convert longs to strings
first_base_offset = _min % 4
last_base_offset = _max % 4
if last_base_offset == 0:
    seq =  seq[first_base_offset:]
else:
    seq =  seq[first_base_offset:(last_base_offset-4)]
    
print "Sequence length:", len(seq)
print "Adjusted sequence:\n", seq.tostring()

Output:

start block: 2 , end block: 25
23 blocks to read.
Sequence in bytes:
AAAAAGCTTTTACTTTTTGCGGCCTAAGTTAGCCAAGCCTAGTAGTTTCTAGAGGCAGAAGTTTTTTTAGTTTCACAGACTCTATTGCGAAG
Sequence length: 90
Adjusted sequence:
AAAGCTTTTACTTTTTGCGGCCTAAGTTAGCCAAGCCTAGTAGTTTCTAGAGGCAGAAGTTTTTTTAGTTTCACAGACTCTATTGCGAAG

Step 4: Apply the blocks (*->N) and masks (upper->lower)

from bisect import bisect_left
# blocks to N       
idx = bisect_left(bstarts, _min) # _min <= self.blockstarts[idx]
while idx < blockcount and bstarts[idx] < _max:
    start = max (bstarts[idx], _min) - _min
    end   = min (bstarts[idx] + bsizes[idx], _max) - _min
    seq[start:end] = array('c', 'N' * (end - start))
    idx +=1 

# masks to lower
idx = bisect_left(mstarts, _min) # _min <= self.maskstarts[idx]
while idx < maskcount and mstarts[idx] < _max:
    start = max (mstarts[idx], _min) - _min
    end   = min (mstarts[idx] + msizes[idx], _max) - _min
    seq[start:end] = array('c', seq[start:end].tostring().lower())
    idx +=1 
print seq.tostring()

Output:

AAAGCTTNNNNNNTTTGCGGCCTAAGTTAGCCAAGCCTAGTAGTTTCTAGAGGCAGAAGTTTTTTNNNNNTCACAGACTCTATTGCGAAG