makewiggle.py

   1 #
   2 #  makewiggle.py
   3 #  ENRAGE
   4 #
   5 import sys
   6 import optparse
   7 import ReadDataset
   8 from commoncode import getConfigParser, getConfigOption, getConfigBoolOption, getConfigIntOption, getConfigFloatOption
   9
  10 print "makewiggle: version 6.8"
  11
  12 try:
  13     import psyco
  14     psyco.full()
  15 except:
  16     print 'psyco not running'
  17
  18
  19 def main(argv=None):
  20     if not argv:
  21         argv = sys.argv
  22
  23     usage = "usage: python %s name rdsfile outfilename [options]"
  24
  25     parser = getParser(usage)
  26     (options, args) = parser.parse_args(argv[1:])
  27
  28     if len(args) < 3:
  29         print usage
  30         sys.exit(1)
  31
  32     name = args[0]
  33     hitfilename = args[1]
  34     outfilename = args[2]
  35
  36     makewiggle(name, hitfilename, outfilename, options.doNormalize, options.color, options.altColor,
  37                options.limitChrom, options.shift, options.doSplit, options.listfilename, options.listPrefix,
  38                options.group, options.startPriority, options.skipRandom, options.withMulti,
  39                options.withSplices, options.doSingle, options.cachePages, options.enforceChr, options.strand,
  40                options.chunk)
  41
  42
  43 def getParser(usage):
  44     parser = optparse.OptionParser(usage=usage)
  45     parser.add_option("--raw", action="store_false", dest="doNormalize")
  46     parser.add_option("--color", dest="color")
  47     parser.add_option("--altcolor", dest="altColor")
  48     parser.add_option("--chrom", dest="limitChrom")
  49     parser.add_option("--shift", type="int", dest="shift")
  50     parser.add_option("--split", action="store_true", dest="doSplit")
  51     parser.add_option("--listfile", dest="listfilename")
  52     parser.add_option("--listprefix", dest="listPrefix")
  53     parser.add_option("--group", dest="group")
  54     parser.add_option("--startPriority", type="float", dest="startPriority")
  55     parser.add_option("--skiprandom", action="store_true", dest="skipRandom")
  56     parser.add_option("--nomulti", action="store_false", dest="withMulti")
  57     parser.add_option("--splices", action="store_true", dest="withSplices")
  58     parser.add_option("--singlebase", action="store_true", dest="doSingle")
  59     parser.add_option("--cache", type="int", dest="cachePages")
  60     parser.add_option("--enforceChr", action="store_true", dest="enforceChr")
  61     parser.add_option("--stranded", dest="strand")
  62     parser.add_option("--maxchunk", type="int", dest="chunk")
  63
  64     configParser = getConfigParser()
  65     section = "makewiggle"
  66     doNormalize = getConfigBoolOption(configParser, section, "doNormalize", True)
  67     color = getConfigOption(configParser, section, "color", None)
  68     altColor = getConfigOption(configParser, section, "altColor", "")
  69     limitChrom = getConfigOption(configParser, section, "limitChrom", None)
  70     shift = getConfigIntOption(configParser, section, "shift", 0)
  71     doSplit = getConfigBoolOption(configParser, section, "doSplit", False)
  72     listfilename = getConfigOption(configParser, section, "listfilename", None)
  73     listPrefix = getConfigOption(configParser, section, "listPrefix", "")
  74     group = getConfigOption(configParser, section, "group", "")
  75     startPriority = getConfigFloatOption(configParser, section, "startPriority", 0.01)
  76     skipRandom = getConfigBoolOption(configParser, section, "skipRandom", False)
  77     withMulti = getConfigBoolOption(configParser, section, "withMulti", True)
  78     withSplices = getConfigBoolOption(configParser, section, "withSplices", False)
  79     doSingle = getConfigBoolOption(configParser, section, "doSingle", False)
  80     cachePages = getConfigIntOption(configParser, section, "cachePages", -1)
  81     enforceChr = getConfigBoolOption(configParser, section, "enforceChr", False)
  82     strand = getConfigOption(configParser, section, "strand", None)
  83     chunk = getConfigIntOption(configParser, section, "chunk", 20)
  84
  85     parser.set_defaults(doNormalize=doNormalize, color=color, altColor=altColor, limitChrom=limitChrom,
  86                         shift=shift, doSplit=doSplit, listfilename=listfilename, listPrefix=listPrefix,
  87                         group=group, startPriority=startPriority, skipRandom=skipRandom, withMulti=withMulti,
  88                         withSplices=withSplices, doSingle=doSingle, cachePages=cachePages, enforceChr=enforceChr,
  89                         strand=strand, chunk=chunk)
  90
  91     return parser
  92
  93
  94 def makewiggle(name, hitfilename, outfilename, doNormalize=True, color=None, altColor="",
  95                limitChrom=None, shift=0, doSplit=False, listfilename=None, listPrefix="",
  96                group="", startPriority=0.01, skipRandom=False, withMulti=True, withSplices=False,
  97                doSingle=False, cachePages=-1, enforceChr=False, strand=None, chunk=20):
  98
  99     priorityIncrement = 0.01
 100     wigType = "bedGraph"
 101
 102     if color is not None:
 103         colorString = " color=%s" % color
 104     else:
 105         colorString = ""
 106
 107     if altColor:
 108         colorString += " altcolor=%s" % altColor
 109
 110     doList = False
 111     if listfilename is not None:
 112         doList = True
 113
 114     chromLimit = False
 115     if limitChrom is not None:
 116         chromLimit = True
 117
 118     if group:
 119         groupName = "group=%s" % group
 120
 121     doCache = False
 122     if cachePages > 0:
 123         doCache = True
 124
 125     maxSpan = chunk * 1000000
 126
 127     isStranded = False
 128     strandedDirection = "both"
 129     if strand is not None:
 130         isStranded = True
 131         if strand == "plus":
 132             strandedDirection = "plusOnly"
 133         elif strand == "minus":
 134             strandedDirection = "minusOnly"
 135
 136         print "will keep track of %s strand(s)" % strandedDirection
 137
 138     if shift:
 139         print "Will shift reads by +/- %d bp according to their sense" % shift
 140         name += "shift=%d" % shift
 141
 142     hitRDS = ReadDataset.ReadDataset(hitfilename, verbose=True, cache=doCache)
 143
 144     if cachePages > hitRDS.getDefaultCacheSize():
 145         hitRDS.setDBcache(cachePages)
 146
 147     readlen = hitRDS.getReadSize()
 148
 149     if doNormalize:
 150         normalizeBy = len(hitRDS) / 1000000.
 151     else:
 152         normalizeBy = 1.
 153
 154     if doList:
 155         listfile = open(listfilename, "w")
 156
 157     priority = startPriority
 158     if not doSplit:
 159         outfile = open(outfilename, "w")
 160         if doList:
 161             listfile.write("%s%s\n" % (listPrefix, outfilename))
 162
 163         outfile.write('track type=%s name="%s" %s priority=%.3f visibility=full%s\n' % (wigType, name, groupName, priority, colorString))
 164
 165     chromList = hitRDS.getChromosomes()
 166     chromList.sort()
 167     for achrom in chromList:
 168         if enforceChr and ("chr" not in achrom):
 169             continue
 170
 171         if chromLimit and achrom != limitChrom:
 172             continue
 173
 174         if skipRandom and "random" in achrom:
 175             continue
 176
 177         if doSplit:
 178             outfile = open("%s.%s" % (outfilename, achrom), "w")
 179             if doList:
 180                 listfile.write("%s%s.%s\n" % (listPrefix, outfilename, achrom))
 181
 182             outfile.write('track type=%s name="%s %s" %s priority=%.3f visibility=full%s\n' % (wigType, name, achrom, groupName, priority, colorString))
 183             priority += priorityIncrement
 184
 185         lastNT = hitRDS.getMaxCoordinate(achrom, doMulti=withMulti, doSplices=withSplices) + readlen
 186         spanStart = 0
 187
 188         previousVal = 0
 189         previousStart = 1
 190         lineIndex = 0
 191         for spanStop in xrange(maxSpan, lastNT+maxSpan, maxSpan):
 192             if spanStop > lastNT:
 193                 spanStop = lastNT
 194
 195             print achrom, spanStart, spanStop
 196             chromModel = hitRDS.getChromProfile(achrom, spanStart, spanStop, withMulti, withSplices, normalizeBy, isStranded, strandedDirection, shiftValue=shift)
 197
 198             for index in xrange(len(chromModel)):
 199                 currentVal = chromModel[index]
 200                 if doSingle:
 201                     outline = "%s %d %.4f\n" % (achrom, spanStart + index, currentVal)
 202                     outfile.write(outline)
 203                     continue
 204
 205                 if currentVal == previousVal:
 206                     continue
 207
 208                 if currentVal != previousVal:
 209                     if previousVal != 0:
 210                         lastpos = index + spanStart
 211                         outline = "%s %d %d %.4f\n" % (achrom, previousStart, lastpos, previousVal)
 212                         outfile.write(outline)
 213                         lineIndex += 1
 214
 215                     previousVal = currentVal
 216                     previousStart = index + spanStart
 217
 218             currentVal = 0
 219             del chromModel
 220             spanStart = spanStop + 1
 221
 222         if doSplit:
 223             outfile.close()
 224
 225         if doSingle:
 226             print index + 1
 227         else:
 228             print lineIndex
 229
 230     if not doSplit:
 231         outfile.close()
 232
 233     if doList:
 234         listfile.close()
 235
 236
 237 if __name__ == "__main__":
 238     main(sys.argv)