10 print 'psyco not running'
14 from commoncode import getMergedRegions, computeRegionBins
16 print '%s: version 2.0' % sys.argv[0]
19 print 'usage: python %s regionfile rdsfile outfilename [-bins numbins] [-field fieldNum] [-raw] [-padregion bp] [-mergeregion bp] [-cache]' % sys.argv[0]
22 regionfilename = sys.argv[1]
24 outfilename = sys.argv[3]
26 if '-raw' in sys.argv:
34 if '-cache' in sys.argv:
38 if '-field' in sys.argv:
39 fieldIndex = sys.argv.index('-field') + 1
40 cField = int(sys.argv[fieldIndex])
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
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
55 if '-bins' in sys.argv:
56 binfield = sys.argv.index('-bins') + 1
57 bins = int(sys.argv[binfield])
59 hitRDS = ReadDataset.ReadDataset(hitfile, verbose=True, cache=doCache)
60 readlen = hitRDS.getReadSize()
61 normalizationFactor = 1.0
63 totalCount = len(hitRDS)
64 normalizationFactor = totalCount / 1000000.
66 chromList = hitRDS.getChromosomes(fullChrom=False)
69 regionDict = getMergedRegions(regionfilename, maxDist=mergeregion, keepLabel=True, verbose=True, chromField=cField, pad=padregion)
71 hitDict = hitRDS.getReadsDict(doMulti=True, findallOptimize=True)
73 (regionsBins, regionsLen) = computeRegionBins(regionsByChromDict, hitDict, bins, readlen, normalizationFactor)
75 outfile = open(outfilename, 'w')
76 for regionID in regionsBins:
78 for binAmount in regionsBins[regionID]:
81 outfile.write('%s\t%s\t%.1f\t%d' % (regionID, regionID, tagCount, Len[gid]))
82 for binAmount in gidBins[gid]:
86 outfile.write('\t%.1f' % (100. * binAmount / tagCount))
88 outfile.write('\t%.1f' % binAmount)