BigWigFile

BigWigFile

• BigWigFile is a class to read BigWig format files.
• It's a python wrapper of Jimmy Kent's C codes. The wrapper code is in external/Kentlib/wWigIO.
• The compiled wWigIO.so should be enough for general use. The BigWigFile class just makes a pythonic interface of the wWigIO.so.

＃ Advantage of this BigWigFile class

1. Easy to understand and convenient to use. Instead of a complicated package, a single wWigIO.so should be enough.
2. Space and time efficient. Speed is very fast. Bigwig files are not opened and closed on each call.
3. Allow to open multiple bigwig files, and multiple handles to one bigwig file.

wWigIO.so

• The wWigIO.so can be found in the installed package, or compiled standalone in the external/Kentlib/wWigIO.
• The wWigIO library can be used directly.

Example: Use wWigIO to open and close bigwig files

import wWigIO

# open two files
wWigIO.open('test.bw')
wWigIO.open('test2.bw')

# open again
wWigIO.open('test.bw')
wWigIO.open('test2.bw')

# close files
wWigIO.close('test.bw')
wWigIO.close('test2.bw')

# close again
wWigIO.close('test.bw')
wWigIO.close('test2.bw'

Output:

open file:test.bw
open file:test2.bw
number of test.bw +1
number of test2.bw +1
number of test.bw -1
number of test2.bw -1
close file(not first node):test.bw
close file(first node):test2.bw

Example: get chrom sizes and fetch wigs by file name

wWigIO.open('test.bw')
chroms = wWigIO.getChromSize('test.bw')
wigs = wWigIO.getIntervals('test.bw', 'chr1', 10, 200)
wWigIO.close('test.bw')
print chrom
print wigs

Output:

{'chr1': 197195432}
[(10, 50, 0.3100000023841858), (50, 100, 0.4099999964237213), (100, 150, 0.33000001311302185), (150, 200, 0.3499999940395355)]

BigWigFile : Make wWigIO pythonic

Although wWigIO.so is enough to use. We'd better to make it more pythonic.

• Make it a Python class with more convenient interface.
• Setup automatic close function to avoid memory leak.

Definition of BigWigFile:

class BigWigFile(object):
    '''
    Fast reader of BigWig format file.
    Usage:
        Open file:
            fh=BigWigFile("test.bw")
        Get chromosome sizes:
            chroms=fh.chromSizes() # {key:value = chrom:size}
        Fetch regions:
            wigs=fh.fetch(chrom="chr1",start=100,stop=200)
            for wig in wigs:
                #do some thing with wig
                print wig
                print wig.chrom,wig.start,wig.stop,wig.score
        Close file:
            fh.close()
        Parameters:
        chrom=None: return empty list.
        start=None: start at first position.
        stop=None:  stop at the end of chromosome.
    '''
   4 lines:     def __init__(self,fname):---------------------------------------------------------
   4 lines:     def chromSizes(self):-------------------------------------------------------------
   9 lines:     def fetch(self,**kwargs):---------------------------------------------------------
  12 lines:     def fetchBed(self,tbed,byPos=False,forcestrand=True):-----------------------------
   3 lines:     def close(self):------------------------------------------------------------------
   3 lines:     def __del__(self): ---------------------------------------------------------------
   4 lines:     def wigToBigWig(wigfile, sizefile, bwfile):---------------------------------------
   4 lines:     def bigWigToWig(bwfile, wigfile):-------------------------------------------------

Example: read bigwig file with BigWigFile class

from ngslib import BigWigFile
fh=BigWigFile("test.bw")
chroms=fh.chromSizes() # {key:value = chrom:size}
print chroms
wigs=fh.fetch(chrom="chr1",start=100,stop=200)
for wig in wigs:
    #do some thing with wig
    print wig
    print wig.chrom,wig.start,wig.stop,wig.score
fh.close() # the bigwig file will be close automatically if we forget to close it.

Output:

open file:test.bw
{'chr1': 197195432}
chr1 100 150 0.330000013113
chr1 99 150 0.330000013113
chr1 150 200 0.34999999404
chr1 149 200 0.34999999404
close file(first node):test.bw

Example: Use Bed as input. Output can be a numpy array if bypos = True. The array can be reversed if forcestrand = True and Bed is on the minus strand.

bwf = BigWigFile("test.bw")
tbed = Bed('chr1\t100\t200')
deptharray = bwf.fetchBed(tbed, bypos=True, forcestrand=False) 
print deptharray[0:10]

Output:

[0.330000013113, 0.330000013113,  ....]

Example: Wiggle <-> BigWig conversion

from ngslib import BigWigFile

# wiggle -> bigwig
BigWigFile.wigToBigWig('test.wig','test.sizes','test.bw') # test.sizes is the chromosome size file with each line like this: "chr1\t197195432"

# bigwig -> wiggle
BigWigFile.bigWigToWig('test.bw','test2.wig')

Then we check the test.wig and test2.wig in shell:

> head test2.wig
variableStep chrom=chr1 span=50
1       0.31
51      0.41
101     0.33
151     0.35
3000901 0.41
3000951 2.34
3001001 2.08
3001051 2.52
3001101 1.41
> head test.wig
variableStep chrom=chr1 span=50
1       0.31
51      0.41
101     0.33
151     0.35
3000901 0.41
3000951 2.34
3001001 2.08
3001051 2.52
3001101 1.41