cistematic/cisstat/analyzego.py

   1 ###########################################################################
   2 #                                                                         #
   3 # C O P Y R I G H T   N O T I C E                                         #
   4 #  Copyright (c) 2003-10 by:                                              #
   5 #    * California Institute of Technology                                 #
   6 #                                                                         #
   7 #    All Rights Reserved.                                                 #
   8 #                                                                         #
   9 # Permission is hereby granted, free of charge, to any person             #
  10 # obtaining a copy of this software and associated documentation files    #
  11 # (the "Software"), to deal in the Software without restriction,          #
  12 # including without limitation the rights to use, copy, modify, merge,    #
  13 # publish, distribute, sublicense, and/or sell copies of the Software,    #
  14 # and to permit persons to whom the Software is furnished to do so,       #
  15 # subject to the following conditions:                                    #
  16 #                                                                         #
  17 # The above copyright notice and this permission notice shall be          #
  18 # included in all copies or substantial portions of the Software.         #
  19 #                                                                         #
  20 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
  21 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
  22 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
  23 # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
  24 # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
  25 # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
  26 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
  27 # SOFTWARE.                                                               #
  28 ###########################################################################
  29 #
  30 # get Gene Ontology Information and distribution
  31
  32 import cistematic.core
  33 from cistematic.core.geneinfo import geneinfoDB
  34 from cistematic.cisstat.score import pvalue
  35 from cistematic.cisstat.helper import hyperGeometric
  36 from math import sqrt
  37 from random import shuffle
  38 import string
  39
  40 def calculateGOStats(geneIDList, fileprefix, excludeIDList = [], roundsRandomization=0, sigLevel = 0.01):
  41     """ calculates GO enrichment (and depletion) statistics for a given list of
  42         geneIDs (assumed to be all from the same genome) using either the
  43         hypergeometric distribution or random sampling for roundsRandomization
  44         rounds if greater than 0. Specific geneID's can be excluded from the
  45         calculations if included in excludeIDList.
  46
  47         Results are saved in files with the given fileprefix. GO Terms that are
  48         larger than 15 genes and that have pvalues < sigLevel / (#num GO categories)
  49         are reported in fileprefix.gosig, whereas genes that matches the GO term are
  50         listed in fileprefix.ZZZZZZZ where the Z's are the GOID.
  51     """
  52     print "calculateGOStats: %s" % fileprefix
  53     if roundsRandomization > 0:
  54         goexpfile = open("%s.goexp" % fileprefix, "w")
  55     gostatfile = open("%s.gostat" % fileprefix, "w")
  56     gozscorefile = open("%s.gozscore" % fileprefix, "w")
  57     gosigfile = open("%s.gosig" % fileprefix, "w")
  58
  59     if len(geneIDList) < 1:
  60         print "Need at the very least one gene!"
  61         return
  62
  63     firstgid = geneIDList[0]
  64     genome = firstgid[0]
  65     cistematic.core.cacheGeneDB(genome)
  66     idb = geneinfoDB(cache=True)
  67     goBin = {}
  68     goPossible = {}
  69     found = {}
  70     locusList = []
  71     excludeLocusList = []
  72     excludeGOBin = {}
  73     zList = []
  74     sigList = []
  75     print "Getting GO list"
  76     goDesc = cistematic.core.allGOTerms(genome)
  77     goInfo = cistematic.core.getAllGOInfo(genome)
  78     print "len(goDesc) = %d" % len(goDesc)
  79     for GOID in goDesc:
  80         goBin[GOID] = 0
  81         goPossible[GOID] = cistematic.core.getGOIDCount(genome, GOID)
  82
  83     for entry in geneIDList:
  84         if entry not in excludeIDList and entry not in locusList:
  85             locusList.append(entry)
  86
  87     for (genome, gID) in excludeIDList:
  88         if gID not in excludeLocusList:
  89             excludeLocusList.append(gID)
  90             try:
  91                 excludeGOTerms = goInfo[(genome, gID)]
  92             except:
  93                 continue
  94
  95             for entry in excludeGOTerms:
  96                 excludeGOTermsFields = entry.split("\t")
  97                 GOID = excludeGOTermsFields[0]
  98                 if GOID not in excludeGOBin:
  99                     excludeGOBin[GOID] = 0
 100
 101                 excludeGOBin[GOID] += 1
 102
 103     print "Getting GO list for locusList"
 104     for locus in locusList:
 105         try:
 106             locusGOTerms = goInfo[locus]
 107         except:
 108             continue
 109
 110         for entry in locusGOTerms:
 111             locusGOTermsFields = entry.split("\t")
 112             GOID = locusGOTermsFields[0]
 113             if GOID not in found:
 114                 found[GOID] = []
 115
 116             if locus not in found[GOID]:
 117                 found[GOID].append(locus)
 118                 goBin[GOID] += 1
 119
 120     numGenes = len(locusList)
 121     goSize = {}
 122     print "Arranging by size"
 123     print "numGenes = %s" % str(numGenes)
 124     for GOID in goBin:
 125         goLen = goBin[GOID]
 126         if goLen not in goSize:
 127             goSize[goLen] = []
 128
 129         goSize[goLen].append(GOID)
 130
 131     goLengths = goSize.keys()
 132     goLengths.sort()
 133     goLengths.reverse()
 134     goList = []
 135     for goLen in goLengths:
 136         for GOID in goSize[goLen]:
 137             goList.append(GOID)
 138
 139     allGIDs = []
 140     try:
 141         theGIDs = cistematic.core.getGenomeGeneIDs(genome)
 142     except:
 143         print "could not get gene entries for %s" % genome
 144
 145     for aGID in theGIDs:
 146         if aGID not in excludeLocusList:
 147             allGIDs.append(aGID)
 148
 149     gokeys = goBin.keys()
 150     mean = {}
 151     standardDev = {}
 152     if roundsRandomization > 0:
 153         print "Get Random sampling"
 154         sample = {}
 155         for sampleNum in range(roundsRandomization):
 156             print "Round %d" % (sampleNum + 1)
 157             sample[sampleNum] = {}
 158             for GOID in gokeys:
 159                 sample[sampleNum][GOID] = 0
 160
 161             shuffle(allGIDs)
 162             sampleGenes = allGIDs[:numGenes]
 163             for gid in sampleGenes:
 164                 try:
 165                     goarray = goInfo[(genome, gid)]
 166                 except:
 167                     continue
 168
 169                 for entry in goarray:
 170                     goarrayFields = entry.split("\t")
 171                     GOID = goarrayFields[0]
 172                     if GOID in gokeys:
 173                         sample[sampleNum][GOID] += 1
 174
 175         print "Calculating stats"
 176         for GOID in gokeys:
 177             mean[GOID] = 0
 178             sumofsquares = 0
 179             for sampleNum in range(roundsRandomization):
 180                 mean[GOID] += sample[sampleNum][GOID]
 181
 182             mean[GOID] = float(mean[GOID]) / float(roundsRandomization)
 183             for sampleNum in range(roundsRandomization):
 184                 sumofsquares += (sample[sampleNum][GOID] - mean[GOID]) * (sample[sampleNum][GOID] - mean[GOID])
 185
 186             standardDev[GOID] = sqrt(sumofsquares / float(roundsRandomization - 1))
 187             goexpfile.write("%s\t%f\t%f\n" % (GOID, mean[GOID], standardDev[GOID]))
 188     else:
 189         # Use hypergeometric
 190         N = float(len(allGIDs))
 191         for GOID in gokeys:
 192             possible = goPossible[GOID]
 193             if GOID in excludeGOBin:
 194                 possible -= excludeGOBin[GOID]
 195
 196             mean[GOID] = (numGenes * possible) / N
 197             standardDev[GOID] = sqrt(numGenes * (possible / N) * (N - possible / N ) * (N - numGenes) / (N - 1.0))
 198
 199     print "Writing out gostat"
 200     for GOID in goList:
 201         count = goBin[GOID]
 202         possible = goPossible[GOID]
 203         if GOID in excludeGOBin:
 204             possible -= excludeGOBin[GOID]
 205
 206         try:
 207             divisor = float(standardDev[GOID])
 208             if divisor > 0:
 209                 zscore = float(count - mean[GOID]) / float(standardDev[GOID])
 210             else:
 211                 zscore = 0.0
 212
 213             if possible > 0:
 214                 percentage = float(count) * 100.0 / float(possible)
 215             else:
 216                 percentage = 0
 217
 218             pval = 1.0
 219             if roundsRandomization > 0:
 220                 if zscore < 0.0:
 221                     pval = pvalue(-1.0 * zscore)
 222                 else:
 223                     pval = pvalue(zscore)
 224             else:
 225                 zscore = count - mean[GOID]
 226                 pval = hyperGeometric(len(allGIDs), possible, numGenes, count)
 227
 228             status = "enriched"
 229             if (count - mean[GOID]) < 0:
 230                 status = "depleted"
 231
 232             gostatfile.write("%s\t%d out of %d\t%2.2f\t%f\t%3.3g\t%s\n" % (GOID, count, possible, percentage, zscore, pval, goDesc[GOID]))
 233             zList.append((zscore, (GOID, count, possible, percentage, zscore, pval, goDesc[GOID])))
 234             if len(goList) > 0 and pval < (sigLevel / float(len(goList))) and possible > 15 and count > 1:
 235                 sigList.append((pval, (GOID, count, possible, percentage, pval, status, goDesc[GOID])))
 236         except:
 237             gostatfile.write("%s\t%d out of %d\t \t%s\n" % (GOID, count, possible, goDesc[GOID]))
 238
 239     gostatfile.close()
 240     print "writing gozscore"
 241     zList.sort()
 242     zList.reverse()
 243     for (zscore, entry) in zList:
 244         gozscorefile.write("%s\t%d out of %d\t%2.2f\t%f\t%3.3g\t%s\n" % entry)
 245
 246     gozscorefile.close()
 247     print "writing gosig"
 248     sigList.sort()
 249     annotCache = {}
 250     for (pval, entry) in sigList:
 251         GOID = entry[0]
 252         gosigfile.write("%s\t%d out of %d\t%2.2f\t%3.3g\t%s\t%s\n" % entry)
 253         goidsigfile = open("%s.%s" % (fileprefix, GOID[3:]), "w")
 254         goidsigfileList = []
 255         if GOID not in found:
 256             print "could not find %s" % GOID
 257             continue
 258
 259         print "%s\t%d genes" % (GOID, len(found[GOID]))
 260         for gID in found[GOID]:
 261             geneSymbol = ""
 262             outsyn = ""
 263             outdesc = ""
 264             geneInfo = ""
 265             if gID in annotCache:
 266                 (geneSymbol, outsyn, outdesc) = annotCache[gID]
 267             else:
 268                 try:
 269                     geneInfo = idb.getGeneInfo(gID)
 270                     geneSymbol = geneInfo[0]
 271                     synonyms = idb.geneIDSynonyms(gID)
 272                     description = idb.getDescription(gID)
 273                     outsyn = string.join(synonyms)
 274                     outdesc = string.join(description, ";")
 275                 except:
 276                     geneSymbol = str(gID)
 277                     try:
 278                         if len(gID) == 2:
 279                             newGeneID = idb.getGeneID(gID[0], gID[1])
 280                             if len(newGeneID) > 1:
 281                                 geneInfo = idb.getGeneInfo(newGeneID)
 282
 283                         if len(geneInfo) > 0:
 284                             geneSymbol = gID[1]
 285                             synonyms = idb.geneIDSynonyms(newGeneID)
 286                             description = idb.getDescription(newGeneID)
 287                             outsyn = string.join(synonyms)
 288                             outdesc = string.join(description, ";")
 289                     except:
 290                         pass
 291
 292             goidsigfileList.append("%s\t%s\t%s\n" % (geneSymbol, outsyn, outdesc))
 293             if gID not in annotCache:
 294                 annotCache[gID] = (geneSymbol, outsyn, outdesc)
 295
 296         goidsigfileList.sort()
 297         for line in goidsigfileList:
 298             goidsigfile.write(line)
 299
 300         goidsigfile.close()
 301
 302     gosigfile.close()
 303     cistematic.core.uncacheGeneDB(genome)