regionBins.py

   1 #
   2 #  regionBins.py
   3 #  ENRAGE
   4 #
   5
   6 try:
   7     import psyco
   8     psyco.full()
   9 except:
  10     print 'psyco not running'
  11
  12 import sys
  13 print '%s: version 2.0' % sys.argv[0]
  14
  15 if len(sys.argv) < 4:
  16     print 'usage: python %s regionfile rdsfile outfilename [-bins numbins] [-field fieldNum] [-raw] [-padregion bp] [-mergeregion bp] [-cache]' % sys.argv[0]
  17     sys.exit(1)
  18
  19 from commoncode import *
  20
  21 regionfilename = sys.argv[1]
  22 hitfile =  sys.argv[2]
  23 outfilename = sys.argv[3]
  24
  25 if '-raw' in sys.argv:
  26     normalize = False
  27     normalizeBins = False
  28 else:
  29     normalize = True
  30     normalizeBins = True
  31
  32 doCache = False
  33 if '-cache' in sys.argv:
  34     doCache = True
  35
  36 cField = 1
  37 if '-field' in sys.argv:
  38     fieldIndex = sys.argv.index('-field') + 1
  39     cField = int(sys.argv[fieldIndex])
  40
  41 padregion = 0
  42 if '-padregion' in sys.argv:
  43     padField = sys.argv.index('-padregion') + 1
  44     padregion = int(sys.argv[padField])
  45     print 'padding %d bp on each side of a region' % padregion
  46
  47 mergeregion = 0
  48 if '-mergeregion' in sys.argv:
  49     mergeField = sys.argv.index('-mergeregion') + 1
  50     mergeregion = int(sys.argv[mergeField])
  51     print 'merging regions closer than %d bp' % mergeregion
  52
  53 bins = 10
  54 if '-bins' in sys.argv:
  55     binfield = sys.argv.index('-bins') + 1
  56     bins = int(sys.argv[binfield])
  57
  58 hitRDS = readDataset(hitfile, verbose = True, cache=doCache)
  59 readlen = hitRDS.getReadSize()
  60 normalizationFactor = 1.0
  61 if normalize:
  62     totalCount = len(hitRDS)
  63     normalizationFactor = totalCount / 1000000.
  64
  65 chromList = hitRDS.getChromosomes(fullChrom=False)
  66 chromList.sort()
  67
  68 regionDict = getMergedRegions(regionfilename, maxDist = mergeregion, keepLabel = True, verbose = True, chromField = cField, pad=padregion)
  69
  70 hitDict = hitRDS.getReadsDict(doMulti=True, findallOptimize=True)
  71
  72 (regionsBins, regionsLen) = computeRegionBins(regionsByChromDict, hitDict, bins, readlen, normalizationFactor)
  73
  74 outfile = open(outfilename, 'w')
  75 for regionID in regionsBins:
  76     tagCount = 0.
  77     for binAmount in regionsBins[regionID]:
  78         tagCount += binAmount
  79     outfile.write('%s\t%s\t%.1f\t%d' % (regionID, regionID, tagCount, Len[gid]))
  80     for binAmount in gidBins[gid]:
  81             if normalizeBins:
  82                 if tagCount == 0:
  83                     tagCount = 1
  84                 outfile.write('\t%.1f' % (100. * binAmount / tagCount))
  85             else:
  86                 outfile.write('\t%.1f' % binAmount)
  87     outfile.write('\n')
  88
  89 outfile.close()