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