makerdsfrombed.py

   1 #
   2 #  makerdsfrombed.py
   3 #  ENRAGE
   4 #
   5 #  Created by Ali Mortazavi on 6/21/08.
   6 #
   7 try:
   8     import psyco
   9     psyco.full()
  10 except:
  11     pass
  12
  13 import sys
  14 import string
  15 import optparse
  16 from commoncode import writeLog, getConfigParser, getConfigIntOption, getConfigBoolOption
  17 import ReadDataset
  18
  19 verstring = "makerdsfrombed: version 2.2"
  20 print verstring
  21
  22
  23 def main(argv=None):
  24     if not argv:
  25         argv = sys.argv
  26
  27     usage = "usage: python %prog label bedfile outrdsfile [--append] [--index] [propertyName::propertyValue] [--cache numPages]"
  28
  29     parser = makeParser(usage)
  30     (options, args) = parser.parse_args(argv[1:])
  31
  32     if len(args) < 3:
  33         print usage
  34         print "\ntreats all imported reads as uniquely mapped\n"
  35         sys.exit(1)
  36
  37     label = args[0]
  38     filename = args[1]
  39     outdbname = args[2]
  40
  41     if options.rnaDataType:
  42         dataType = "RNA"
  43     else:
  44         dataType = "DNA"
  45
  46     propertyList = []
  47     for arg in args:
  48         if "::" in arg:
  49             (pname, pvalue) = arg.strip().split("::")
  50             propertyList.append((pname, pvalue))
  51
  52     makerdsfrombed(label, filename, outdbname, options.init, dataType, options.doIndex, options.cachePages, propertyList)
  53
  54
  55 def makeParser(usage=""):
  56     parser = optparse.OptionParser(usage=usage)
  57     parser.add_option("--append", action="store_false", dest="init")
  58     parser.add_option("--index", action="store_true", dest="doIndex")
  59     parser.add_option("--cache", type="int", dest="cachePages")
  60     parser.add_option("--RNA", action="store_true", dest="rnaDataType")
  61
  62     configParser = getConfigParser()
  63     section = "makerdsfrombed"
  64     init = getConfigBoolOption(configParser, section, "init", True)
  65     rnaDataType = getConfigBoolOption(configParser, section, "RNA", False)
  66     doIndex = getConfigBoolOption(configParser, section, "doIndex", False)
  67     cachePages = getConfigIntOption(configParser, section, "cachePages", 100000)
  68
  69     parser.set_defaults(init=init, rnaDataType=rnaDataType, doIndex=doIndex, cachePages=cachePages)
  70
  71     return parser
  72
  73
  74 def makerdsfrombed(label, filename, outdbname, init=True, dataType="DNA", doIndex=False, cachePages=100000, propertyList=[]):
  75     readsize = 0
  76     padsize = 0
  77     index = 0
  78     insertSize = 100000
  79
  80     writeLog(outdbname + ".log", verstring, string.join(sys.argv[1:]))
  81
  82     infile = open(filename,"r")
  83
  84     rds = ReadDataset.ReadDataset(outdbname, init, dataType, verbose=True)
  85     if not init:
  86         rds.dropIndex()
  87
  88     #check that our cacheSize is better than the dataset's default cache size
  89     defaultCacheSize = rds.getDefaultCacheSize()
  90     if cachePages > defaultCacheSize:
  91         if init:
  92             rds.setDBcache(cachePages, default=True)
  93         else:
  94             rds.setDBcache(cachePages)
  95
  96     if len(propertyList) > 0:
  97         rds.insertMetadata(propertyList)
  98
  99     insertList = []
 100     for line in infile:
 101         if "track" in line:
 102             continue
 103
 104         fields = line.split()
 105         if readsize == 0:
 106             readsize = abs(int(fields[1]) - int(fields[2]))
 107             if init:
 108                 rds.insertMetadata([("readsize", readsize+1)])
 109                 rds.insertMetadata([("imported_from_bed", "True")])
 110
 111         chrom = fields[0]
 112         start = int(fields[1])
 113         stop = int(fields[2])
 114         sense = fields[5]
 115         readID = "%s-%s" % (label, str(index))
 116         insertList.append((readID, chrom, start, stop, sense, 1.0, "", ""))
 117         if index % insertSize == 0:
 118             rds.insertUniqs(insertList)
 119             insertList = []
 120             print ".",
 121             sys.stdout.flush()
 122
 123         index += 1
 124
 125     if len(insertList) > 0:
 126         rds.insertUniqs(insertList)
 127
 128     countString = "%d unique reads" % index
 129     print countString
 130
 131     writeLog(outdbname + ".log", verstring, countString)
 132
 133     if doIndex:
 134         print "building index...."
 135         if cachePages > defaultCacheSize:
 136             rds.setDBcache(cachePages)
 137             rds.buildIndex(cachePages)
 138         else:
 139             rds.buildIndex(defaultCacheSize)
 140
 141
 142 if __name__ == "__main__":
 143     main(sys.argv)