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", action="store_false", 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     tabfile = open(tabFileName)
  56     previous = ""
  57
  58     posList = []
  59     for line in tabfile:
  60         fields = line.split("\t")
  61         current = fields[0]
  62         if previous == current:
  63             continue
  64
  65         previous = current
  66         chrom = fields[1][3:]
  67         posList.append((chrom, (int(fields[2]) + int(fields[3]))/2))
  68
  69     feats = featuresIntersecting("human", posList, radius, cistype)
  70     featkeys = feats.keys()
  71     featkeys.sort()
  72     for (chrom, pos) in featkeys:
  73         print "chr%s:%d-%d\t%s" % (chrom, pos, pos + 20, str(feats[(chrom, pos)]))
  74
  75
  76 if __name__ == "__main__":
  77     main(sys.argv)