#!/usr/bin/env python2.7importsysimportbx.bbi.bigwig_fileimportnumpyasnpimportmatplotlibmatplotlib.use("Agg")importmatplotlib.pyplotaspltimportosimportstructdefgetChromosomeSizesFromBigWig(bwname):csize={}fh=open(os.path.expanduser(bwname),"rb")magic=fh.read(4)ifmagic=='&\xfc\x8f\x88':endianness='<'elifmagic=='\x88\x8f\xfc&':endianness='>'else:raiseIOError("The file is not in bigwig format")(version,zoomLevels,chromosomeTreeOffset,fullDataOffset,fullIndexOffset,fieldCount,definedFieldCount,autoSqlOffset,totalSummaryOffset,uncompressBufSize,reserved)=struct.unpack(endianness+'HHQQQHHQQIQ',fh.read(60))ifversion<3:raiseIOError("Bigwig files version <3 are not supported")fh.seek(chromosomeTreeOffset)magic=fh.read(4)ifmagic=='\x91\x8c\xcax':endianness='<'elifmagic=='x\xca\x8c\x91':endianness='>'else:raiseValueError("Wrong magic for this bigwig data file")(blockSize,keySize,valSize,itemCount,reserved)=struct.unpack(endianness+'IIIQQ',fh.read(28))(isLeaf,reserved,count)=struct.unpack(endianness+'BBH',fh.read(4))forninrange(count):(key,chromId,chromSize)=struct.unpack(endianness+str(keySize)+'sII',fh.read(keySize+2*4))csize[key.replace('\x00','')]=chromSizereturncsize# These functions were copied from simpleHilbertCurve (https://github.com/dentearl/simpleHilbertCurve.git)defd2xy(n,d):""" take a d value in [0, n**2 - 1] and map it to an x, y value (e.g. c, r). """assert(d<=n**2-1)t=dx=y=0s=1while(s<n):rx=1&(t/2)ry=1&(t^rx)x,y=rot(s,x,y,rx,ry)x+=s*rxy+=s*ryt/=4s*=2returnx,ydefrot(n,x,y,rx,ry):""" rotate/flip a quadrant appropriately """ifry==0:ifrx==1:x=n-1-xy=n-1-yreturny,xreturnx,ySample1=bx.bbi.bigwig_file.BigWigFile(open(sys.argv[1],'rb'))Sample2=bx.bbi.bigwig_file.BigWigFile(open(sys.argv[2],'rb'))csize=getChromosomeSizesFromBigWig(sys.argv[1])imsize=640n=imsize-1forchromincsize.keys():sys.stderr.write("Processing chromosome %s\n"%chrom)chrom_length=csize[chrom]sys.stderr.write("Processing Sample1\n")hc=np.zeros((imsize,imsize))data=Sample1.get_as_array(chrom,0,chrom_length)data[np.isnan(data)]=0foriinxrange(0,len(chrom_length)):c,r=d2xy(n,int(round((n**2-1)*i/chrom_length)))hc[r][c]+=data[i]# dump hilbert curve in image and arraynp.save("Sample1_%s"%chrom,hc)plt.imshow(hc,origin=0,alpha=0.5,cmap=plt.cm.Blues)sys.stderr.write("Processing Sample2\n")hc=np.zeros((imsize,imsize))data=Sample2.get_as_array(chrom,0,chrom_length)data[np.isnan(data)]=0foriinxrange(0,len(chrom_length)):c,r=d2xy(n,int(round((n**2-1)*i/chrom_length)))hc[r][c]+=data[i]# dump hilbert curve in image and arraynp.save("Sample2_%s"%chrom,hc)plt.imshow(hc,origin=0,alpha=0.5,cmap=plt.cm.Reds)plt.title("Hilbert curve for %s"%chrom)plt.savefig("HC_%s.png"%chrom)plt.close()