geneStartBins.py

   1 #
   2 #  geneStartBins.py
   3 #  ENRAGE
   4 #
   5
   6 try:
   7     import psyco
   8     psyco.full()
   9 except:
  10     pass
  11
  12 # originally from version 1.3 of geneDownstreamBins.py
  13 from commoncode import *
  14 from cistematic.genomes import Genome
  15 from cistematic.core.geneinfo import geneinfoDB
  16 import sys
  17
  18 print '%s: version 2.0' % sys.argv[0]
  19 if len(sys.argv) < 4:
  20     print 'usage: python %s genome rdsfile outfilename [-max regionSize] [-raw] [-cache]' % sys.argv[0]
  21     print '\n\twhere regionSize is the optional maximum region in bp\n'
  22     sys.exit(1)
  23
  24 genome = sys.argv[1]
  25 hitfile =  sys.argv[2]
  26 outfilename = sys.argv[3]
  27
  28 standardMinDist = 3000
  29 if '-max' in sys.argv:
  30     standardMinDist = int(sys.argv[sys.argv.index('-max') + 1])
  31
  32 if '-raw' in sys.argv:
  33     normalize = False
  34     normalizeBins = False
  35 else:
  36     normalize = True
  37     normalizeBins = True
  38
  39 doCache = False
  40 if '-cache' in sys.argv:
  41     doCache = True
  42
  43 bins = 10
  44 standardMinThresh = standardMinDist / bins
  45
  46 hitRDS = readDataset(hitfile, verbose = True, cache=doCache)
  47 readlen = hitRDS.getReadSize()
  48 normalizationFactor = 1.0
  49 if normalize:
  50     totalCount = len(hitRDS)
  51     normalizationFactor = totalCount / 1000000.
  52
  53 hg = Genome(genome)
  54 idb = geneinfoDB(cache=True)
  55
  56 gidDict = {}
  57 geneinfoDict = idb.getallGeneInfo(genome)
  58 featuresDict = hg.getallGeneFeatures()
  59
  60 #infile = open(infilename)
  61 outfile = open(outfilename,'w')
  62
  63 gidList = hg.allGIDs()
  64 gidList.sort()
  65 for gid in gidList:
  66     symbol = 'LOC' + gid
  67     geneinfo = ''
  68     featureList = []
  69     try:
  70         geneinfo = geneinfoDict[gid]
  71         featureList = featuresDict[gid]
  72         symbol = geneinfo[0][0]
  73     except:
  74         print geneinfo
  75     newfeatureList = []
  76     if len(featureList) == 0:
  77         continue
  78     for (ftype, chrom, start, stop, fsense) in featureList:
  79         if (start, stop) not in newfeatureList:
  80             newfeatureList.append((start, stop))
  81     if chrom not in hitDict:
  82         continue
  83     newfeatureList.sort()
  84     if len(newfeatureList) < 1:
  85         #print '%s %s %d' % (gid, symbol, -1)
  86         #outfile.write('%s\t%s\t%d\n' % (gid, symbol, -1))
  87         continue
  88     glen = standardMinDist / 2
  89     if fsense == 'F':
  90         nextGene = hg.leftGeneDistance((genome, gid), glen * 2)
  91         if nextGene < glen * 2:
  92                 glen = nextGene / 2
  93         if glen < 1:
  94                 glen = 1
  95         gstart = newfeatureList[0][0] - glen
  96         if gstart < 0:
  97                 gstart = 0
  98         gstop = newfeatureList[0][0]  + glen
  99     else:
 100         nextGene = hg.rightGeneDistance((genome, gid), glen * 2)
 101         if nextGene < glen * 2:
 102             glen = nextGene / 2
 103         if glen < 1:
 104             glen = 1
 105         gstart = newfeatureList[-1][1] - glen
 106         gstop = newfeatureList[-1][1] + glen
 107     tagCount = 0
 108     if glen < standardMinDist / 2:
 109         continue
 110     binList = [0] * bins
 111     for (tagStart, sense, weight) in hitDict[chrom]:
 112         tagStart -= gstart
 113         if tagStart >= 2 * glen:
 114             break
 115         if tagStart > 0:
 116             tagCount += weight
 117             if fsense == 'R':
 118                 # we are relying on python's integer division quirk
 119                 binID = tagStart / standardMinThresh
 120                 binList[binID] += weight
 121             else:
 122                 rdist = 2 * glen - tagStart
 123                 binID = rdist / standardMinThresh
 124                 binList[binID] += weight
 125     if tagCount < 2:
 126         continue
 127     print '%s %s %d %d %s' % (gid, symbol, tagCount, glen, str(binList))
 128     outfile.write('%s\t%s\t%d\t%d' % (gid, symbol, tagCount, glen))
 129     for binAmount in binList:
 130         outfile.write('\t%d' % binAmount)
 131     outfile.write('\n')
 132 #infile.close()
 133 outfile.close()
 134