featureIntersects.py

   1 #
   2 #  featureIntersects.py
   3 #  ENRAGE
   4 #
   5
   6 try:
   7     import psyco
   8     psyco.full()
   9 except:
  10     pass
  11
  12 import sys
  13 import optparse
  14 from cistematic.core import featuresIntersecting
  15 from commoncode import getConfigParser, getConfigOption, getConfigIntOption
  16
  17
  18 print "featureIntersects: version 1.1"
  19
  20
  21 def main(argv=None):
  22     if not argv:
  23         argv = sys.argv
  24
  25     usage = "usage: python %s tabfile [--cistype type] [--radius radius]"
  26
  27     parser = makeParser(usage)
  28     (options, args) = parser.parse_args(argv[1:])
  29
  30     if len(args) < 1:
  31         print usage
  32         sys.exit(1)
  33
  34     tabfile = args[0]
  35
  36     featureIntersects(tabfile, options.cistype, options.radius)
  37
  38
  39 def makeParser(usage=""):
  40     parser = optparse.OptionParser(usage=usage)
  41     parser.add_option("--cistype", dest="cistype")
  42     parser.add_option("--radius", type="int", dest="radius")
  43
  44     configParser = getConfigParser()
  45     section = "featureIntersects"
  46     cistype = getConfigOption(configParser, section, "cistype", "TFBSCONSSITES")
  47     radius = getConfigIntOption(configParser, section, "radius", 100)
  48
  49     parser.set_defaults(cistype=cistype, radius=radius)
  50
  51     return parser
  52
  53
  54 def featureIntersects(tabFileName, cistype="TFBSCONSSITES", radius=100):
  55
  56     posList = getPositionList(tabFileName)
  57     feats = featuresIntersecting("hsapiens", posList, radius, cistype)
  58     featkeys = feats.keys()
  59     featkeys.sort()
  60     for (chrom, pos) in featkeys:
  61         print "chr%s:%d-%d\t%s" % (chrom, pos, pos + 20, str(feats[(chrom, pos)]))
  62
  63
  64 def getPositionList(tabFileName):
  65     previous = ""
  66     posList = []
  67     tabfile = open(tabFileName)
  68     for line in tabfile:
  69         fields = line.split("\t")
  70         current = fields[0]
  71         if previous == current:
  72             continue
  73
  74         previous = current
  75         chrom = fields[1][3:]
  76         posList.append((chrom, (int(fields[2]) + int(fields[3]))/2))
  77
  78     return posList
  79
  80
  81 if __name__ == "__main__":
  82     main(sys.argv)