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