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     groupName = ""
 119     if group:
 120         groupName = "group=%s" % group
 121
 122     doCache = False
 123     if cachePages > 0:
 124         doCache = True
 125
 126     maxSpan = chunk * 1000000
 127
 128     isStranded = False
 129     strandedDirection = "both"
 130     if strand is not None:
 131         isStranded = True
 132         if strand == "plus":
 133             strandedDirection = "plusOnly"
 134         elif strand == "minus":
 135             strandedDirection = "minusOnly"
 136
 137         print "will keep track of %s strand(s)" % strandedDirection
 138
 139     if shift:
 140         print "Will shift reads by +/- %d bp according to their sense" % shift
 141         name += "shift=%d" % shift
 142
 143     hitRDS = ReadDataset.ReadDataset(hitfilename, verbose=True, cache=doCache)
 144
 145     if cachePages > hitRDS.getDefaultCacheSize():
 146         hitRDS.setDBcache(cachePages)
 147
 148     readlen = hitRDS.getReadSize()
 149
 150     if doNormalize:
 151         normalizeBy = len(hitRDS) / 1000000.
 152     else:
 153         normalizeBy = 1.
 154
 155     if doList:
 156         listfile = open(listfilename, "w")
 157
 158     priority = startPriority
 159     if not doSplit:
 160         outfile = open(outfilename, "w")
 161         if doList:
 162             listfile.write("%s%s\n" % (listPrefix, outfilename))
 163
 164         outfile.write('track type=%s name="%s" %s priority=%.3f visibility=full%s\n' % (wigType, name, groupName, priority, colorString))
 165
 166     chromList = hitRDS.getChromosomes()
 167     chromList.sort()
 168     for achrom in chromList:
 169         if enforceChr and ("chr" not in achrom):
 170             continue
 171
 172         if chromLimit and achrom != limitChrom:
 173             continue
 174
 175         if skipRandom and "random" in achrom:
 176             continue
 177
 178         if doSplit:
 179             outfile = open("%s.%s" % (outfilename, achrom), "w")
 180             if doList:
 181                 listfile.write("%s%s.%s\n" % (listPrefix, outfilename, achrom))
 182
 183             outfile.write('track type=%s name="%s %s" %s priority=%.3f visibility=full%s\n' % (wigType, name, achrom, groupName, priority, colorString))
 184             priority += priorityIncrement
 185
 186         lastNT = hitRDS.getMaxCoordinate(achrom, doMulti=withMulti, doSplices=withSplices) + readlen
 187         spanStart = 0
 188
 189         previousVal = 0
 190         previousStart = 1
 191         lineIndex = 0
 192         for spanStop in xrange(maxSpan, lastNT+maxSpan, maxSpan):
 193             if spanStop > lastNT:
 194                 spanStop = lastNT
 195
 196             print achrom, spanStart, spanStop
 197             chromModel = hitRDS.getChromProfile(achrom, spanStart, spanStop, withMulti, withSplices, normalizeBy, isStranded, strandedDirection, shiftValue=shift)
 198
 199             for index in xrange(len(chromModel)):
 200                 currentVal = chromModel[index]
 201                 if doSingle:
 202                     outline = "%s %d %.4f\n" % (achrom, spanStart + index, currentVal)
 203                     outfile.write(outline)
 204                     continue
 205
 206                 if currentVal == previousVal:
 207                     continue
 208
 209                 if currentVal != previousVal:
 210                     if previousVal != 0:
 211                         lastpos = index + spanStart
 212                         outline = "%s %d %d %.4f\n" % (achrom, previousStart, lastpos, previousVal)
 213                         outfile.write(outline)
 214                         lineIndex += 1
 215
 216                     previousVal = currentVal
 217                     previousStart = index + spanStart
 218
 219             currentVal = 0
 220             del chromModel
 221             spanStart = spanStop + 1
 222
 223         if doSplit:
 224             outfile.close()
 225
 226         if doSingle:
 227             print index + 1
 228         else:
 229             print lineIndex
 230
 231     if not doSplit:
 232         outfile.close()
 233
 234     if doList:
 235         listfile.close()
 236
 237
 238 if __name__ == "__main__":
 239     main(sys.argv)