1. Davide Cittaro
  2. utils

Source

utils / plot_map.py

import numpy as np
import os
import struct
import bx.bbi.bigwig_file
import scipy.ndimage
import scipy.interpolate
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import sys


def getChromosomeSizesFromBigWig(bwname):
  csize = {}
  fh = open(os.path.expanduser(bwname), "rb")
  magic = fh.read(4)
  if magic == '&\xfc\x8f\x88':
    endianness = '<'
  elif magic == '\x88\x8f\xfc&':
    endianness = '>'
  else:
    raise IOError("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))
  if version < 3:
    raise IOError("Bigwig files version <3 are not supported")
  fh.seek(chromosomeTreeOffset)
  magic = fh.read(4)
  if magic == '\x91\x8c\xcax':
    endianness = '<'
  elif magic == 'x\xca\x8c\x91':
    endianness = '>'
  else:
    raise ValueError("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))
  for n in range(count):
    (key, chromId, chromSize) = struct.unpack(endianness + str(keySize) + 'sII', fh.read(keySize + 2 * 4))
    csize[key.replace('\x00', '')] = chromSize
  return csize


Sample1 = 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 = 640

image = np.zeros((imsize,imsize))

for chrom in csize.keys():
  if '_' in chrom: continue
  if chrom == 'chrX' or chrom =='chrY' or chrom =='chrM': continue
  print(chrom)
  ar1 = Sample1.get_as_array(chrom, 0, csize[chrom])
  ar1[np.isnan(ar1)] = 0
  ar1[ar1 > imsize] = imsize
  ar2 = Sample2.get_as_array(chrom, 0, csize[chrom])
  ar2[np.isnan(ar2)] = 0
  ar2[ar2 > imsize] = imsize
  xbin = np.linspace(np.min(ar1), np.max(ar1), 641)
  ybin = np.linspace(np.min(ar2), np.max(ar2), 641)
  t_image, xd, yd = np.histogram2d(ar1,ar2, bins=[xbin, ybin])
  image += t_image

np.save(sys.argv[3], image)
plt.imshow(np.log(image), origin=0)
plt.yticks(np.linspace(0, 640, 3), np.linspace(np.min(ar1),np.max(ar1), 3))
plt.xticks(np.linspace(0, 640, 3), np.linspace(np.min(ar2),np.max(ar2), 3))
plt.xlabel(sys.argv[2])
plt.ylabel(sys.argv[1])
plt.title(sys.argv[3])
plt.savefig(sys.argv[3] + '.png')