getgosig.py

   1 try:
   2     import psyco
   3     psyco.full()
   4 except:
   5     pass
   6
   7 from cistematic.genomes import Genome
   8 from math import log
   9 import os.path
  10 import sys
  11 import optparse
  12 import matplotlib
  13 from pylab import *
  14 from commoncode import getConfigParser, getConfigIntOption
  15
  16
  17 print "getgosig: version 2.2"
  18
  19 def main(argv=None):
  20     if not argv:
  21         argv = sys.argv
  22
  23     usage = "usage: python %prog genome outimage gofileroot1 title1 cohortsize1 [gofileroot2 title2 cohortsize2 ...] [--fontsize pts] [--length in] [--width in]"
  24
  25     parser = makeParser(usage)
  26     (options, args) = parser.parse_args(argv[1:])
  27
  28     if len(args) < 5:
  29         print usage
  30         sys.exit(1)
  31
  32     genome = args[0]
  33     imagename =  args[1]
  34
  35     conditionList = args[2:]
  36     conditions = len(conditionList) / 3
  37     fileroots = []
  38     titles = []
  39     for index in range(conditions):
  40         conditionIndex = index * 3
  41         fileroots.append(conditionList[conditionIndex])
  42         titles.append((conditionList[conditionIndex + 1], "(%s)" % conditionList[conditionIndex + 2]))
  43
  44     getgosig(genome, imagename, fileroots, titles, options.fontSize, options.length, options.width)
  45
  46
  47 def makeParser(usage=""):
  48     parser = optparse.OptionParser(usage=usage)
  49     parser.add_option("--fontsize", type="int", dest="fontSize")
  50     parser.add_option("--length", type="int", dest="length")
  51     parser.add_option("--width", type="int", dest="width")
  52
  53     configParser = getConfigParser()
  54     section = "getgosig"
  55     fontSize = getConfigIntOption(configParser, section, "fontSize", 5)
  56     length = getConfigIntOption(configParser, section, "length", 10)
  57     width = getConfigIntOption(configParser, section, "width", 7)
  58
  59     parser.set_defaults(fontSize=fontSize, length=length, width=width)
  60
  61     return parser
  62
  63
  64 def getgosig(genome, imagename, fileroots=[], titles=[], fontSize=5, length=10, width=7):
  65     hg = Genome(genome)
  66     allgodesc = hg.allGOterms()
  67     godesc = []
  68
  69     matplotlib.use("Agg")
  70
  71     doGray = False
  72
  73     rootdir = "./"
  74     htmlname = imagename[:-4] + ".html"
  75
  76     ceiling = 40.0
  77     goterms = []
  78     goscores = {}
  79     numgenes = {}
  80     possiblegenes = {}
  81     flatArray = []
  82
  83     highestPval = 0.0
  84     lowestPval = 1.0
  85     for sigfile in fileroots:
  86         infile = open(rootdir + sigfile + ".gosig", "r")
  87         for line in infile:
  88             if "depleted" in line:
  89                 continue
  90
  91             fields = line.split("\t")
  92             if fields[0] not in goterms:
  93                 goterms.append(fields[0])
  94                 goscores[fields[0]] = []
  95                 numgenes[fields[0]] = []
  96                 possiblegenes[fields[0]] = 0
  97
  98             if float(fields[3]) > highestPval:
  99                 highestPval = float(fields[3])
 100
 101             if float(fields[3]) < lowestPval:
 102                 lowestPval = float(fields[3])
 103
 104     print highestPval
 105     print lowestPval
 106
 107     boundaryScore = score = -1 * log(highestPval) /  (2.0 * ceiling) + 0.49
 108     print boundaryScore
 109
 110     cdict = {"red": ((0.0, 1.0, 1.0), (boundaryScore, 1.0, 0.1), (1.0, 1.0, 1.0)),
 111              "green": ((0.0, 1.0, 1.0), (boundaryScore, 1.0, 0.1), (1.0, 1.0, 1.0)),
 112              "blue": ((0.0, 1.0, 1.0), (boundaryScore, 1.0, 0.75), (1.0, 0.0, 0.0))
 113     }
 114
 115     mymap = matplotlib.colors.LinearSegmentedColormap("my_colormap", cdict, 1024)
 116
 117     goindex = 0
 118     for zfile in fileroots:
 119         infile = open(rootdir + zfile + ".gozscore", "r")
 120         for line in infile:
 121             fields = line.split()
 122             goindex += 1
 123             if fields[0] not in goterms:
 124                 continue
 125
 126             score = -1 * log(float(fields[7])) /  (2.0 * ceiling)
 127             if score < -0.5:
 128                 score = -0.5
 129
 130             if score > 0.5:
 131                 score = 0.5
 132
 133             score += 0.5
 134             if doGray:
 135                 score = 1 - score
 136
 137             goscores[fields[0]].append(score)
 138             numgenes[fields[0]].append(fields[1])
 139             possiblegenes[fields[0]] = int(fields[4])
 140
 141     goindex /= len(fileroots)
 142
 143     gokeys = goscores.keys()
 144     gosortarray = []
 145     for term in gokeys:
 146         gosortarray.append(goscores[term] + [term])
 147
 148     gosortarray.sort()
 149
 150     htmlfile = open(htmlname, "w")
 151     htmlfile.write('<html><head><title>GO Analysis</title></head><body><table border="1">')
 152     htmlfile.write("<tr><th>Description</th><th>possible</th>")
 153     for entry in titles:
 154         htmlfile.write("<th>%s<br>%s</th>" % entry)
 155
 156     htmlfile.write("</tr>\n")
 157     tableLines = []
 158
 159     for entry in gosortarray:
 160         term = entry[-1]
 161         outline = "%s:\t" % term
 162         for entry in goscores[term]:
 163             outline += str(round(entry, 4)) + "\t"
 164
 165         print outline
 166         htmlLine = "<tr><th>%s</th><th>%d</th>" % (allgodesc[term], possiblegenes[term])
 167         index = 0
 168         for fileroot in fileroots:
 169             gofile = fileroot + "." + term[3:]
 170             ngene = numgenes[term][index]
 171             if os.path.exists(gofile):
 172                 htmlLine += '<td><a href="%s">%s</a></td>' % (gofile, ngene)
 173             else:
 174                 htmlLine += "<td>%s</td>" % (ngene)
 175
 176             index += 1
 177
 178         tableLines.append(htmlLine + "</tr>\n")
 179         flatArray.append(goscores[term])
 180         godesc.append(allgodesc[term])
 181
 182     tableLines.reverse()
 183     for line in tableLines:
 184         htmlfile.write(line)
 185
 186     htmlfile.write("<tr><th>Cohort Size:</th>")
 187     htmlfile.write("</tr>\n")
 188     htmlfile.write("</table></body></html>")
 189
 190     figure(figsize=(length, width))
 191     myaxe = axes([0.3, 0.1, 0.55, 0.75])
 192
 193     Z = array(flatArray)
 194     print Z.shape
 195     if doGray:
 196         c = pcolor(Z, cmap=cm.gray, vmin=0.0, vmax=1.0)
 197     else:
 198         c = pcolor(Z, cmap=mymap, vmin=0.0, vmax=1.0)
 199
 200     c.set_linewidth(0.1)
 201     clim(0.0, 1.0)
 202
 203     ind = arange(len(fileroots))
 204     width = 0.5
 205
 206     coordy = 0.1
 207     deltaX = 1.0
 208     deltaY = 1.0
 209
 210     pcolorAxes = c.get_axes()
 211     for entry in gosortarray:
 212         term = entry[-1]
 213         coordx = 0.4
 214         for genenum in numgenes[term]:
 215             if len(genenum) == 1:
 216                 genenum = "    " + genenum
 217             elif len(genenum) == 2:
 218                 genenum = "  " + genenum
 219
 220             pcolorAxes.text(coordx, coordy, genenum, fontsize=fontSize)
 221             coordx += deltaX
 222
 223         coordy += deltaY
 224
 225     coordx = 0
 226     for (line1,line2) in titles:
 227         pcolorAxes.text(coordx + 0.1, coordy + 3 * deltaY + 0.5, line1, fontsize=int(fontSize*1.5))
 228         pcolorAxes.text(coordx + 0.1, coordy + deltaY, line2, fontsize=int(fontSize*1.5))
 229         coordx += deltaX
 230
 231     setp(gca(), "xticks", [])
 232     setp(gca(), "xticklabels", [])
 233     setp(gca(), "yticks", arange(len(godesc)))
 234     setp(gca(), "yticklabels", godesc)
 235     locs, labels = yticks()
 236     setp(labels, fontsize=fontSize)
 237     setp(labels, verticalalignment="bottom")
 238     setp(gca(), "ylim", [0, len(godesc)])
 239
 240     figtext(0.3,0.02, str(goindex - len(gokeys)) + " additional GO Terms below threshold of significance", fontsize=fontSize*2)
 241
 242     d = colorbar(orientation="vertical", drawedges=False)
 243     for t in d.ax.get_yticklabels():
 244         t.set_fontsize(0)
 245
 246     locs, labels = yticks()
 247     setp(labels, fontsize=5)
 248     pcolorAxes.text(conditions + 1,len(godesc), str(lowestPval), fontsize=fontSize*2)
 249     pcolorAxes.text(conditions + 1,boundaryScore * len(godesc), str(highestPval), fontsize=fontSize*2)
 250
 251     savefig(imagename, dpi=250)
 252     show()
 253
 254
 255 if __name__ == "__main__":
 256     main(sys.argv)