cistematic/experiments/analyzeMotifs.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 # basic analysis code
  31 from cistematic.core.motif import matrixRow
  32 import cistematic.core
  33 import string
  34
  35 try:
  36     from pysqlite2 import dbapi2 as sqlite
  37 except:
  38     from sqlite3 import dbapi2 as sqlite
  39 from os import environ
  40
  41 if environ.get("CISTEMATIC_ROOT"):
  42     cisRoot = environ.get("CISTEMATIC_ROOT")
  43 else:
  44     cisRoot = "/proj/genome"
  45
  46 knownMotifsDB = "%s/db/known_motifs.db" % cisRoot
  47
  48
  49 class AnalyzeMotifs:
  50     motifSize = {}
  51     coappear  = {}
  52     motifToAnnotations = {}
  53     annotationEntries = []
  54     genomeMatches = {}
  55     loadedGenomes = []
  56     dbName = ""
  57     analysisID = "default"
  58
  59
  60     def initializeKnownMotifs(self, source=""):
  61         """  loads known motifs from motif database, optionally limiting by source.
  62         """
  63         self.annotationEntries = []
  64         dbkm = sqlite.connect(knownMotifsDB, timeout=60)
  65         sqlkm = dbkm.cursor()
  66         if len(source) > 0:
  67             sqlkm.execute("select * from motifs where source = '%s' " % source)
  68         else:
  69             sqlkm.execute("select * from motifs")
  70
  71         dbEntries = sqlkm.fetchall()
  72         dbkm.close()
  73
  74         for entry in dbEntries:
  75             (index, source, name, theseq, reference, other_info) = entry
  76             self.annotationEntries.append((index, source, name, theseq, reference, other_info))
  77
  78
  79     def findAnnotations(self, motID, thresh=-1, numberN=0):
  80         mot = self.findMotif(motID)
  81         if thresh < 0:
  82             thresh = self.threshold
  83
  84         for entry in self.annotationEntries:
  85             (id, source, name, theseq, reference, other) = entry
  86             loc = mot.locateMotif(theseq, thresh, numberN)
  87             if len(loc) > 0 :
  88                 self.motifToAnnotations[motID].append((source, name, theseq, reference, other, loc))
  89
  90             if len(self.motifToAnnotations[motID]) >= 500:
  91                 break
  92
  93
  94     def findAnnotationsRE(self, motID, numMismatches=0):
  95         """ Calls the motif's consensus locator on each entry in annotationEntries.
  96         """
  97         mot = self.findMotif(motID)
  98         if numMismatches > 0:
  99             mot.initializeMismatchRE(numMismatches)
 100         else:
 101             mot.initializeRE()
 102
 103         for entry in self.annotationEntries:
 104             (id, source, name, annotSeq, reference, other) = entry
 105             # Must use non-RE optimized version for short annotations
 106             if len(annotSeq) < len(mot):
 107                 pad = "N" * (len(mot) - len(annotSeq))
 108                 annotSeq = pad + annotSeq + pad
 109                 loc = mot.locateConsensus(annotSeq)
 110             else:
 111                 loc = mot.locateConsensusRE(annotSeq)
 112
 113             if len(loc) > 0 :
 114                 self.motifToAnnotations[motID].append((source, name, annotSeq, reference, other, loc))
 115
 116             if len(self.motifToAnnotations[motID]) >= 500:
 117                 break
 118
 119
 120     def printAnnotations(self):
 121         """ prints the motif annotations for every motif result.
 122         """
 123         for mot in self.getResults():
 124             annotations =   self.motifToAnnotations[mot.tagID]
 125             print "motif %s\t%s" % (mot.tagID, mot.buildConsensus())
 126             for annotation in annotations:
 127                 (source, name, annotSeq, reference, other, loc) = annotation
 128                 print "%s:%s\t%s\t%s\t%s" % (source, name, reference, annotSeq, str(loc))
 129                 print other
 130             print
 131
 132
 133     def printAnnotationsShort(self):
 134         """ prints a compressed version of the annotations.
 135         """
 136         for motID in self.motifToAnnotations.keys():
 137             for annotation in self.motifToAnnotations[motID]:
 138                 print "%s: %s (%s)\t%s - %s" % (annotation[0], annotation[1], annotation[4], annotation[2], annotation[5])
 139
 140
 141     def returnAnnotations(self, motID):
 142         """ returns the [annotations] for a particular motID
 143         """
 144         try:
 145             return self.motifToAnnotations[motID]
 146         except KeyError:
 147             return []
 148
 149
 150     def annotateMotifs(self, thresh=0.0, numberN=0):
 151         self.mlog( "Annotating Motifs with threshold %d and %d extra Ns" % (1.0 + thresh, numberN))
 152         if len(self.annotationEntries) == 0:
 153             self.initializeKnownMotifs()
 154
 155         for mot in self.getResults():
 156             mot.setThreshold(thresh)
 157             self.motifToAnnotations[mot.tagID] = []
 158             self.findAnnotations(mot.tagID, thresh, numberN)
 159
 160
 161     def annotateConsensus(self, numMismatches=0, source=""):
 162         self.mlog( "Annotating Consensus")
 163         if len(self.annotationEntries) == 0:
 164             self.initializeKnownMotifs(source)
 165
 166         for mot in self.getResults():
 167             self.motifToAnnotations[mot.tagID] = []
 168             self.findAnnotationsRE(mot.tagID, numMismatches)
 169
 170
 171     def printConsensus(self):
 172         """ print the consensus for every motif result.
 173         """
 174         for mot in self.getResults():
 175             print "motif %s\t%s" % (mot.tagID, mot.buildConsensus())
 176
 177
 178     def formatPWM(self, aPWM):
 179         """ format a PWM into a printable string.
 180         """
 181         aRow = ""
 182         cRow = ""
 183         gRow = ""
 184         tRow = ""
 185         for col in aPWM:
 186             aRow = string.join([aRow, str(round(col[matrixRow["A"]],4))], "\t")
 187             cRow = string.join([cRow, str(round(col[matrixRow["C"]],4))], "\t")
 188             gRow = string.join([gRow, str(round(col[matrixRow["G"]],4))], "\t")
 189             tRow = string.join([tRow, str(round(col[matrixRow["T"]],4))], "\t")
 190
 191         formattedPWM =  "A:\t%s\nC:\t%s\nG:\t%s\nT:\t%s\n" % (aRow, cRow, gRow, tRow)
 192
 193         return formattedPWM
 194
 195
 196     def appendGeneToMotif(self,mTagID, geneID, pos):
 197         """ make an entry in the geneToMotif table.
 198         """
 199         (genome, locus) = geneID
 200         (loc, sense) = pos
 201         stmt = "INSERT into geneToMotif(ID, expID, analysisID, mTagID, genome, locus, location, sense) values (NULL, '%s', '%s',  '%s', '%s', '%s', '%s', '%s')" % (self.experimentID, self.analysisID, mTagID, genome, locus, loc, sense)
 202         res = self.sql(stmt, "commit")
 203
 204
 205     def appendGeneToMotifBatch(self, batch):
 206         """ make a batch of entries in the geneToMotif table.
 207         """
 208         batchInserts = []
 209         stmt = "INSERT into geneToMotif(ID, expID, analysisID, mTagID, genome, locus, location, sense) values (NULL, ?, ?,  ?, ?, ?, ?, ?)"
 210         for entry in batch:
 211             (mTagID, geneID, pos) = entry
 212             (genome, locus) = geneID
 213             (loc, sense) = pos
 214             batchInserts.append((self.experimentID, self.analysisID, mTagID, genome, locus, loc, sense))
 215
 216         res = self.batchsql(stmt, batchInserts)
 217
 218
 219     def geneToMotifKeys(self, thekey="geneID"):
 220         """ return the keys to the geneToMotif table. The default key is geneID, otherwise returns mTagID.
 221         """
 222         results = []
 223         if thekey == "geneID":
 224             stmt = "SELECT distinct genome, locus from geneToMotif where expID = '%s' and analysisID = '%s' order by locus" % (self.experimentID, self.analysisID)
 225         else:
 226             stmt = "SELECT distinct mTagID from geneToMotif where expID = '%s' and analysisID = '%s' order by mTagID" % (self.experimentID, self.analysisID)
 227
 228         res = self.sql(stmt)
 229         for entry in res:
 230             if thekey == "geneID":
 231                 (genome, locus) = entry
 232                 results.append((str(genome), str(locus)))
 233             else:
 234                 mTagID = entry[0]
 235                 results.append(str(mTagID))
 236
 237         return results
 238
 239
 240     def appendMotifNeighbors(self,mTagID, match, condition="", geneEntry=""):
 241         """ make an entry in the motifNeighbors table.
 242         """
 243         (chromo, pos) = match
 244         (genome, chromNum) = chromo
 245         (loc, mSense) = pos
 246         if geneEntry != "":
 247             (start, stop, sense, geneID) = geneEntry
 248             (genome2, locus) = geneID
 249         else:
 250             start = "-"
 251             stop = "-"
 252             sense = "-"
 253             locus = "NO GENE IN RADIUS"
 254
 255         stmt = "INSERT into motifNeighbors(ID, expID, analysisID, mTagID, genome, chromNum, location, motifSense, start, stop, sense, locus, condition) values (NULL, '%s','%s',  '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s')" % (self.experimentID, self.analysisID, mTagID, genome, chromNum, loc, mSense, start, stop, sense, locus, condition)
 256         res = self.sql(stmt, "commit")
 257
 258
 259     def appendMotifNeighborsBatch(self, batch):
 260         """ make a batch of entries in the motifNeighbors table.
 261         """
 262         batchInserts = []
 263         stmt = "INSERT into motifNeighbors(ID, expID, analysisID, mTagID, genome, chromNum, location, motifSense, start, stop, sense, locus, condition) values (NULL, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)"
 264         for entry in batch:
 265             (mTagID, match, condition, geneEntry) = entry
 266             (chromo, pos) = match
 267             (genome, chromNum) = chromo
 268             (loc, mSense) = pos
 269             if geneEntry != "":
 270                 (start, stop, sense, geneID) = geneEntry
 271                 (genome2, locus) = geneID
 272             else:
 273                 start = "-"
 274                 stop = "-"
 275                 sense = "-"
 276                 locus = "NO GENE IN RADIUS"
 277
 278             batchInserts.append((self.experimentID, self.analysisID, mTagID, genome, chromNum, loc, mSense, start, stop, sense, locus, condition))
 279
 280         res = self.batchsql(stmt, batchInserts)
 281
 282
 283     def motifNeighborsKeys(self, condition=""):
 284         """ returns a [list of motif ID's] in the motifNeighbor table. Can restrict to a particular condition.
 285         """
 286         results = []
 287         if condition != "":
 288             stmt = "SELECT distinct mTagID from motifNeighbors where expID = '%s' and analysisID = '%s' and condition ='%s' order by mTagID" % (self.experimentID, self.analysisID, condition)
 289         else:
 290             stmt = "SELECT distinct mTagID from motifNeighbors where expID = '%s' and analysisID = '%s' order by mTagID" % (self.experimentID, self.analysisID)
 291
 292         res = self.sql(stmt)
 293         for entry in res:
 294             mTagID = entry[0]
 295             results.append(mTagID)
 296
 297         return results
 298
 299
 300     def motifNeighbors(self,mTagID, condition="", discardNullMatches=False):
 301         """ get entries in the geneToMotif table that match a particular Motif & condition.
 302         """
 303         results = []
 304         genome = ""
 305         chromNum = ""
 306         loc = ""
 307         mSense = ""
 308         start = ""
 309         stop = ""
 310         sense = ""
 311         locus = ""
 312         stmt = "select distinct genome, chromNum, location, motifSense, start, stop, sense, locus, condition from motifNeighbors where expID='%s' and analysisID = '%s' and mTagID='%s' " % (self.experimentID, self.analysisID, mTagID)
 313         if discardNullMatches:
 314             stmt += " and condition <> 'NONE' "
 315
 316         if condition != '':
 317             stmt +=  " and condition = '%s' " % (condition)
 318
 319         stmt += " order by ID desc "
 320         res = self.sql(stmt)
 321
 322         for entry in res:
 323             (genome, chromNum, loc, mSense, start, stop, sense, locus, condition) = entry
 324             match = ((genome, chromNum),(int(loc), mSense))
 325             geneEntry = (start, stop, sense, (genome, locus))
 326             results.append((match, geneEntry, condition))
 327
 328         return results
 329
 330
 331     def motifToGene(self, mTagID):
 332         """ returns a list of matching geneIDs with location and sense found for a given motif.
 333         """
 334         results = []
 335         stmt = "SELECT distinct genome, locus, location, sense from geneToMotif where expID = '%s' and analysisID = '%s' and mTagID = '%s' order by locus" % (self.experimentID, self.analysisID, mTagID)
 336         res = self.sql(stmt)
 337         for entry in res:
 338             (genome, locus, location, sense) = entry
 339             results.append(((genome, locus), (int(location), sense)))
 340
 341         return results
 342
 343
 344     def geneToMotif(self, geneID):
 345         """ returns a list of matching motifs with location and sense found for a given geneID
 346         """
 347         results = []
 348         (genome, locus) = geneID
 349         stmt = "SELECT distinct mTagID, location, sense from geneToMotif where expID = '%s' and analysisID = '%s' and genome = '%s' and locus = '%s' order by location" % (self.experimentID, self.analysisID, genome, locus)
 350         res = self.sql(stmt)
 351         for entry in res:
 352             (mTagID, location, sense) = entry
 353             results.append((mTagID, (location, sense)))
 354
 355         return results
 356
 357
 358     def mapMotifs(self, Threshold=90.0, numberN=0, runList=[], ignoreList=[], enforceSanity=True, mapAllSeqs=True, verbose=True):
 359         """ find occurences of result motifs in the current genepool using PWMs. Slow.
 360         """
 361         currentGenome = ""
 362         currentDataset = -1
 363         genepoolKeys = self.genepool.keys()
 364         posResults = []
 365         for mot in self.getResults():
 366             mTagID = mot.tagID
 367             runID = mTagID.split("-")[0]
 368             (pName, dataID, setID, tStamp, motArray) = self.getRun(int(runID))
 369             if enforceSanity and not mot.isSane():
 370                 self.mlog("mapMotifs: not mapping %s - failed sanity check" % (mTagID))
 371                 continue
 372
 373             if mot.tagID in ignoreList:
 374                 self.mlog("mapMotifs: not mapping %s" % mot.tagID)
 375                 continue
 376
 377             if len(runList) > 0:
 378                 if int(runID) not in runList:
 379                     self.mlog("mapMotifs: not mapping run %s" % runID)
 380                     continue
 381
 382             if mapAllSeqs:
 383                 dataID = 0
 384
 385             self.mlog("mapMotifs: mapping %s with threshold %f and at most %d N" % (mTagID, Threshold, numberN))
 386             if dataID <> currentDataset:
 387                 currentDataset = dataID
 388                 for geneID in self.getSetting(dataID):
 389                     if geneID[0] <> currentGenome:
 390                         currentGenome = geneID[0]
 391
 392                     try:
 393                         if geneID not in genepoolKeys:
 394                             self.genepool[geneID] = cistematic.core.retrieveSeq(geneID, self.upstream, self.cds, self.downstreamself.geneDB, False, self.boundToNextGene)
 395                             genepoolKeys.append[geneID]
 396                     except:
 397                         self.mlog("mapMotifs: could not load %s" % (str(geneID)))
 398                         break
 399
 400             for geneID in genepoolKeys:
 401                 if verbose:
 402                     print "mapMotifs: %s\n" % str(geneID)
 403
 404                 motifPos = mot.locateMotif(self.genepool[geneID], Threshold, numberN)
 405
 406                 if len(motifPos) > 0:
 407                     for pos in motifPos:
 408                         posResults.append((mTagID, geneID, pos))
 409
 410         self.appendGeneToMotifBatch(posResults)
 411
 412
 413     def mapConsensus(self, runList=[], ignoreList=[], enforceSanity=True, mapAllSeqs=True, numMismatches=0):
 414         """ find occurences of result motifs in the current genepool using regular expressions, allowing
 415             for a number of mismatches.
 416         """
 417         currentGenome = ""
 418         currentDataset = -1
 419         genepoolKeys = self.genepool.keys()
 420         posResults = []
 421         for mot in self.getResults():
 422             mTagID = mot.tagID
 423             runID = mTagID.split("-")[0]
 424             (pName, dataID, setID, tStamp, motArray) = self.getRun(int(runID))
 425             if mot.tagID in ignoreList:
 426                 self.mlog("mapConsensus: not mapping %s" % mot.tagID)
 427                 continue
 428
 429             if len(runList) > 0:
 430                 if int(runID) not in runList:
 431                     self.mlog("mapConsensus: not mapping run %s" % runID)
 432                     continue
 433
 434             if enforceSanity and not mot.isSane():
 435                 self.mlog("mapConsensus: not mapping %s - failed sanity check" % (mTagID))
 436                 continue
 437
 438             if mapAllSeqs:
 439                 dataID = 0
 440
 441             if numMismatches > 0:
 442                 mot.initializeMismatchRE(numMismatches)
 443             else:
 444                 mot.initializeRE()
 445             self.mlog("mapConsensus: mapping %s with %s mismatches" % (mTagID, numMismatches))
 446             if dataID <> currentDataset:
 447                 currentDataset = dataID
 448                 for geneID in self.getSetting(dataID):
 449                     if geneID[0] <> currentGenome:
 450                         currentGenome = geneID[0]
 451
 452                     try:
 453                         if geneID not in genepoolKeys:
 454                             self.genepool[geneID] = cistematic.core.retrieveSeq(geneID, self.upstream, self.cds, self.downstream)
 455                             genepoolKeys.append[geneID]
 456                     except:
 457                         self.mlog("mapConsensus: could not load %s" % (str(geneID)))
 458                         break
 459
 460             for geneID in genepoolKeys:
 461                 motifPos = mot.locateConsensusRE(self.genepool[geneID])
 462
 463                 if len(motifPos) > 0:
 464                     for pos in motifPos:
 465                         posResults.append((mTagID, geneID, pos))
 466
 467         self.appendGeneToMotifBatch(posResults)
 468
 469
 470     def mapFeatures(self, radius):
 471         """ returns features within a certain radius in bp of all matches.
 472         """
 473         FeatureList = []
 474         for mot in self.getResults():
 475             mTagID = mot.tagID
 476             self.mlog("mapFeatures: mapping %s using a radius of %d bp" % (mTagID, radius))
 477             for match in self.motifToGene(mTagID):
 478                 matchFeatures = cistematic.core.retrieveFeatures(match, radius)
 479                 for entry in matchFeatures:
 480                     FeatureList.append((mTagID, match, entry))
 481
 482         return FeatureList
 483
 484
 485     def limitGeneEntries(self, geneList, lowerBound, upperBound):
 486         results = []
 487         for entry in geneList:
 488             if entry[1] < lowerBound or entry[0] > upperBound:
 489                 continue
 490
 491             results.append(entry)
 492
 493         return results
 494
 495
 496     def mapNeighbors(self, radius, annotate=True):
 497         """ returns genes on a chromosome within a certain radius in bp.
 498         """
 499         localGeneList = {}
 500         prevChromo = ""
 501         neighborResults = []
 502         for mot in self.getResults():
 503             mTagID = mot.tagID
 504             motLen = len(mot)
 505             motRadius = motLen / 2
 506             mtgList = self.motifToGene(mTagID)
 507             self.mlog("mapNeighbors: mapping %s using a radius of %d bp (%d instances)" % (mTagID, radius, len(mtgList)))
 508             index = 0
 509             for match in mtgList:
 510                 (chromo, hit) = match
 511                 if annotate and chromo != prevChromo:
 512                     prevChromo = chromo
 513                     localGeneList = cistematic.core.getChromoGeneEntries(chromo)
 514
 515                 index += 1
 516                 if (index % 1000) == 0:
 517                     print "."
 518
 519                 matchCounter = 0
 520                 matchCDS = []
 521                 if annotate:
 522                     (chromo, hit) = match
 523                     matchpos = int(hit[0])
 524                     lowerBound = matchpos - radius
 525                     upperBound = matchpos + radius
 526                     match2 = (chromo, (matchpos + motRadius, hit[1]))
 527                     matchFeatures = cistematic.core.retrieveFeatures(match2, motRadius, "CDS")
 528                     for entry in matchFeatures:
 529                         matchCDS.append((chromo[0], entry[0]))
 530
 531                     geneEntriesList = self.limitGeneEntries(localGeneList, lowerBound, upperBound)
 532                     for geneEntry in geneEntriesList:
 533                         beg = int(geneEntry[0])
 534                         end = int(geneEntry[1])
 535                         sense = geneEntry[2]
 536                         gID = geneEntry[3]
 537                         if gID in matchCDS: # matching within the coding sequence
 538                             neighborResults.append((mTagID, match, "CDS", geneEntry))
 539                             matchCounter += 1
 540                         elif matchpos >= beg and matchpos <= end: # not in CDS, but in gene
 541                             neighborResults.append((mTagID, match, "GENE", geneEntry))
 542                             matchCounter += 1
 543                         elif matchpos < beg:
 544                             if sense == "F":
 545                                 neighborResults.append((mTagID, match, "UPSTREAM", geneEntry))
 546                             else:
 547                                 neighborResults.append((mTagID, match, "DOWNSTREAM", geneEntry))
 548
 549                             matchCounter += 1
 550                         else:
 551                             if sense == "F":
 552                                 neighborResults.append((mTagID, match, "DOWNSTREAM", geneEntry))
 553                             else:
 554                                 neighborResults.append((mTagID,match, "UPSTREAM", geneEntry))
 555
 556                             matchCounter += 1
 557
 558                 if matchCounter < 1:
 559                     neighborResults.append((mTagID, match, "NONE", ""))
 560
 561         self.appendMotifNeighborsBatch(neighborResults)
 562
 563
 564     def printMotifToGene(self, motifList=[], runList=[]):
 565         motKeys = self.geneToMotifKeys(thekey="mTagID")
 566         if len(motifList) == 0 and len(runList) == 0:
 567             for tag in motKeys:
 568                 print "Motif %s is found in: %s" % (tag, str(self.motifToGene(tag)))
 569         else:
 570             for tag in motKeys:
 571                 runID = tag.split("-")[0]
 572                 if len(runList) > 0:
 573                     if runID not in runList:
 574                         continue
 575
 576                 if tag in motifList:
 577                     print "Motif %s is found in: %s" % (tag, str(self.motifToGene(tag)))
 578
 579
 580     def motifToGeneStat(self):
 581         tags = self.geneToMotifKeys(thekey="mTagID")
 582         counter = 0
 583         min = len(self.motifToGene(tags[0]))
 584         max = 0
 585         for tag in tags:
 586             numGenes = len(self.motifToGene(tag))
 587             print "%s: %d" % (tag, numGenes)
 588             if numGenes > max:
 589                 max = numGenes
 590             elif numGenes < min:
 591                 min = numGenes
 592
 593             counter += numGenes
 594
 595         print "for a total of %d matches - min: %d max: %d" % (counter, min, max)
 596
 597
 598     def printGeneToMotif(self, geneList = []):
 599         if len(geneList) == 0:
 600             for geneID in self.geneToMotifKeys():
 601                 print "Gene %s has the following motifs: %s" % (str(geneID), str(self.geneToMotif(geneID)))
 602         else:
 603             for geneID in self.geneToMotifKeys():
 604                 if geneID in geneList:
 605                     print "Gene %s has the following motifs: %s" % (str(geneID), str(self.geneToMotif(geneID)))
 606
 607
 608     def geneToMotifStat(self):
 609         genes = self.geneToMotifKeys()
 610         counter = 0
 611         min = len(self.geneToMotif(genes[0]))
 612         max = 0
 613         for gene in genes:
 614             numMotifs = len(self.geneToMotif(gene))
 615             print "%s - %s: %d" % (gene[0], gene[1], numMotifs)
 616             if numMotifs > max:
 617                 max = numMotifs
 618             elif numMotifs < min:
 619                 min = numMotifs
 620
 621             counter += numMotifs
 622
 623         print "for a total of %d matches - min: %d max: %d" % (counter, min, max)
 624
 625
 626     def printGeneProfile(self, geneID):
 627         print "\nMotif matches for %s " % (str(geneID))
 628         geneProfile = self.makeGeneProfile(geneID)
 629         positions = geneProfile.keys()
 630         if len(positions) > 0:
 631             positions.sort()
 632             for pos in positions:
 633                 print "%s %s" % (str(pos), str(geneProfile[pos]))
 634         else:
 635             print "Gene had no matching motifs"
 636
 637
 638     def makeGeneProfile(self, geneID):
 639         geneBucket = {}
 640         geneMotifs = self.geneToMotif(geneID)
 641         if len(geneMotifs) > 0:
 642             for mot in geneMotifs:
 643                 (motID, pos) = mot
 644                 (loc, sense) = pos
 645                 if int(loc) not in geneBucket.keys():
 646                     geneBucket[int(loc)] = []
 647
 648                 geneBucket[int(loc)].append(mot)
 649
 650         return geneBucket
 651
 652
 653     def getMotifMatches(self, motifIDs=[]):
 654         results = {}
 655         if len(motifIDs) < 1:
 656             motifIDs = self.geneToMotifKeys(thekey="mTagID")
 657
 658         for motID in motifIDs:
 659             results[motID] = []
 660             mot = self.findMotif(motID)
 661             motLength = len(mot)
 662             for match in self.motifToGene(motID):
 663                 (chrom, hit) = match
 664                 (pos, sense) = hit
 665                 if sense == "F":
 666                     results[motID].append(self.genepool[chrom][pos:pos + motLength])
 667                 else:
 668                     results[motID].append(cistematic.core.complement(self.genepool[chrom][pos:pos + motLength], motLength))
 669
 670         return results
 671
 672
 673     def summarizeMotifs(self, fileName):
 674         """ saves the number of occurences and actual occurence PWMs of motifs to fileName.
 675         """
 676         motDict = {}
 677         motText = []
 678         try:
 679             if len(fileName) < 1:
 680                 raise IOError
 681
 682             outFile = open(fileName, "a")
 683             self.mlog("Saving motif summary")
 684             for motID in self.geneToMotifKeys(thekey="mTagID"):
 685                 matchNum = 0
 686                 mot = self.findMotif(motID)
 687                 motLength = len(mot)
 688                 motDict[motID] = []
 689                 for index in range(motLength):
 690                     motDict[motID].append([0.0, 0.0, 0.0, 0.0])
 691
 692                 for match in self.motifToGene(motID):
 693                     matchNum += 1
 694                     (chromo, hit) = match
 695                     (pos, sense) = hit
 696                     if sense == "F":
 697                         matchSeq = self.genepool[chromo][pos:pos + motLength]
 698                     else:
 699                         matchSeq = cistematic.core.complement(self.genepool[chromo][pos: pos + motLength], motLength)
 700
 701                     for index in range(motLength):
 702                         NT = matchSeq[index]
 703                         NT = NT.upper()
 704                         if NT in ["A", "C", "G", "T"]:
 705                             motDict[motID][index][matrixRow[NT]] += 1
 706
 707                 motLine = "motif %s\t %s matches\t'%s'\n %s\n" % (motID, str(matchNum), mot.buildConsensus(), self.formatPWM(motDict[motID]))
 708                 motText.append(motLine)
 709
 710             outFile.writelines(motText)
 711             outFile.close()
 712         except:
 713             self.mlog("Could not save motif summary to file %s\n" % fileName)
 714
 715
 716     def printMotifNeighbors(self):
 717         for motID in self.motifNeighborsKeys():
 718             print "Matches for motif %s" % motID
 719             currentHit = ()
 720             for entry in self.motifNeighbors(motID):
 721                 if entry[0] != currentHit:
 722                     print "====================="
 723                     print "Match %s" % str(entry[0])
 724                     currentHit = entry[0]
 725
 726                 print "\tGene %s:" % str(entry[1])
 727                 try:
 728                     goInfo = cistematic.core.getGOInfo(entry[1][3])
 729                     for entry in goInfo:
 730                         print "\t\t%s" % str(entry)
 731                 except:
 732                     pass
 733
 734                 print "----------------"
 735
 736
 737     def summarizeMotifNeighbors(self, fileName):
 738         """ saves the number of occurences and PWMs of motifs within the gene
 739             neighborhood radius as mapped using mapNeighbors() to fileName.
 740         """
 741         motDict = {}
 742         motText = []
 743         try:
 744             if len(fileName) < 1:
 745                 raise IOError
 746             outFile = open(fileName, "a")
 747             self.mlog("Saving neighbor summary")
 748             for motID in self.motifNeighborsKeys():
 749                 matchNum = 0
 750                 mot = self.findMotif(motID)
 751                 motLength = len(mot)
 752                 motDict[motID] = []
 753                 for index in range(motLength):
 754                     motDict[motID].append([0.0, 0.0, 0.0, 0.0])
 755
 756                 currentMatch = ""
 757                 for entry in self.motifNeighbors(motID, discardNullMatches = True):
 758                     if currentMatch != entry[0]:
 759                         currentMatch = entry[0]
 760                         (genechrom, loc) = entry[0]
 761                         (pos, sense) = loc
 762                         geneID = entry[1][3]
 763                         (geno, gene) = geneID
 764                         if gene != "NO GENE IN RADIUS":
 765                             matchNum += 1
 766                             if sense == "F":
 767                                 matchSeq = self.genepool[genechrom][pos:pos + motLength]
 768                             else:
 769                                 matchSeq = cistematic.core.complement(self.genepool[genechrom][pos: pos + motLength], motLength)
 770
 771                             for index in range(motLength):
 772                                 NT = matchSeq[index]
 773                                 NT = NT.upper()
 774                                 if NT in ["A", "C", "G", "T"]:
 775                                     motDict[motID][index][matrixRow[NT]] += 1
 776
 777                 motLine = "motif %s\t %s matches\n %s\n" % (motID, str(matchNum), self.formatPWM(motDict[motID]))
 778                 motText.append(motLine)
 779
 780             outFile.writelines(motText)
 781             outFile.close()
 782         except:
 783             self.mlog("Could not save motif neighbors to file %s\n" % fileName)
 784
 785
 786     def saveMotifNeighbors(self, fileName, fullAnnotations=True):
 787         """ save every occurence of the motifs with any adjoining gene(s) to fileName.
 788             Records annotations and GO terms if available when fullAnnotations=True.
 789         """
 790         goDict = {}
 791         annotDict = {}
 792         currentGenome = ""
 793         neighborList = []
 794         self.mlog("Saving motif neighbors to file %s\n" % fileName)
 795         if True:
 796             if len(fileName) < 1:
 797                 raise IOError
 798             outFile = open(fileName, "a")
 799             for motID in self.motifNeighborsKeys():
 800                 matchNum = 0
 801                 mot = self.findMotif(motID)
 802                 motLength = len(mot)
 803                 currentMatch = ""
 804                 for entry in self.motifNeighbors(motID):
 805                     neighborList = []
 806                     if currentMatch != entry[0]:
 807                         matchNum += 1
 808                         currentMatch = entry[0]
 809                         (genechrom, loc) = entry[0]
 810                         (genome, chromo) = genechrom
 811                         (pos, sense) = loc
 812                         if fullAnnotations and genome != currentGenome:
 813                             currentGenome = genome
 814                             goDict = cistematic.core.getAllGOInfo(genome)
 815                             annotDict = cistematic.core.getAllAnnotInfo(genome)
 816
 817                     start = entry[1][0]
 818                     stop = entry[1][1]
 819                     geneSense = entry[1][2]
 820                     geneID = entry[1][3]
 821                     (geno, gene) = geneID
 822                     condition = entry[2]
 823                     if sense == "F":
 824                         matchSeq = self.genepool[genechrom][pos:pos + motLength]
 825                     else:
 826                         matchSeq = cistematic.core.complement(self.genepool[genechrom][pos: pos + motLength], motLength)
 827
 828                     currentEntry = "%s\t%i\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s" % (motID, matchNum, genome, chromo, pos, sense, condition, start, stop, geneSense, geno, gene, matchSeq)
 829                     if fullAnnotations:
 830                         try:
 831                             goInfo = goDict.get(geneID, [])
 832                             annotInfo = annotDict.get(geneID,[])
 833                             annotDescription = ""
 834                             prevGoInfo = ""
 835                             if self.debugSQL:
 836                                 print "%s\t%s\t%s" % (str(geneID), str(goInfo), str(annotInfo))
 837                             for annot in annotInfo:
 838                                 annotDescription += annot + " "
 839                             if len(goInfo) > 0:
 840                                 for entry in goInfo:
 841                                     if entry == prevGoInfo:
 842                                         continue
 843                                     neighborList.append("%s\t%s\t%s\n" % (currentEntry, annotDescription, entry))
 844                                     prevGoInfo = entry
 845                             else:
 846                                 neighborList.append("%s\t%s\n" % (currentEntry, annotDescription))
 847                         except:
 848                             neighborList.append("%s\n" % currentEntry)
 849                     else:
 850                         neighborList.append("%s\n" % currentEntry)
 851
 852                     outFile.writelines(neighborList)
 853
 854             outFile.close()
 855         else:
 856             self.mlog("Could not save motif neighbors to file %s\n" % fileName)
 857
 858         return neighborList
 859
 860
 861     def buildMotifSize(self):
 862         mots = self.getResults()
 863         index = 0
 864         for mot in mots:
 865             motLength = len(mot)
 866             if motLength not in self.motifSize.keys():
 867                 self.motifSize[motLength] =[]
 868
 869             self.motifSize[motLength].append((index, mot.buildConsensus()))
 870             index += 1
 871
 872
 873     def findIdenticalMotifs(self):
 874         # FIXME!!!!
 875         identicalMotifs = []
 876
 877         for motLength in self.motifSize.keys():
 878             motList = self.motifSize[motLength]
 879             motListLen = len(motList)
 880             print "doing motif length %d" % motLength
 881             for pos in range(motListLen):
 882                 index = pos + 1
 883                 while index < motListLen:
 884                     if motList[pos][1] == motList[index][1]:
 885                         identicalMotifs.append((self.results[pos].tagID, self.results[index].tagID))
 886
 887                     index += 1
 888
 889         return identicalMotifs
 890
 891
 892     def createAnalysis(self, dbName=""):
 893         """ creates the analysis SQL tables. Should only be run once per database file.
 894         """
 895         if len(dbName) == 0:
 896             dbName = self.expFile
 897
 898         db = sqlite.connect(dbName, timeout=60)
 899         try:
 900             sql = db.cursor()
 901             sql.execute("CREATE table geneToMotif(ID INTEGER PRIMARY KEY, expID varchar, analysisID varchar, mTagID varchar, genome varchar, locus varchar, location varchar, sense varchar)")
 902             sql.execute("CREATE table genomeMatches(ID INTEGER PRIMARY KEY, expID varchar, analysisID varchar, mTagID varchar, genome varchar, locus varchar, location varchar, sense varchar, threshold varchar)")
 903             sql.execute("CREATE table motifNeighbors(ID INTEGER PRIMARY KEY, expID varchar, analysisID varchar, mTagID varchar, genome varchar, chromNum varchar, location varchar, motifSense varchar, start varchar, stop varchar, sense varchar, locus varchar, condition varchar)")
 904             sql.execute("CREATE index mot1 on geneToMotif(expID, analysisID, mTagID)")
 905             sql.execute("CREATE index mot2 on genomeMatches(expID, analysisID, mTagID)")
 906             sql.execute("CREATE index mot3 on motifNeighbors(expID, analysisID, mTagID)")
 907             sql.execute("CREATE index locus1 on geneToMotif(locus)")
 908             sql.execute("CREATE index locus2 on motifNeighbors(locus)")
 909             sql.execute("CREATE index condition1 on motifNeighbors(condition)")
 910             db.commit()
 911             sql.close()
 912             db.close()
 913             self.mlog("Created analysis tables in database %s" % dbName)
 914         except:
 915             db.close()
 916             self.mlog("Could not create tables in database %s" % dbName)
 917             self.mlog("WARNING: perhaps you have called createAnalysis twice on the same database?")
 918
 919
 920     def loadAnalysis(self, analysisID="default", dbName=""):
 921         """ changes the analysisID to use the specified one, or use default. Must be used before reading or
 922             writing any data to analysis tables.
 923         """
 924         self.analysisID = analysisID
 925         if len (dbName) > 0:
 926             self.dbName = dbName
 927         else:
 928             self.dbName = self.expFile
 929
 930
 931     def resetAnalysis(self):
 932         """ currently zeros out some of the analysis structures. obsolete ?
 933         """
 934         self.mlog("resetting analysis")
 935         self.motifSize = {}
 936         self.motifToAnnotations = {}
 937         self.coappear  = {}
 938
 939
 940     def saveAnalysis(self):
 941         """ currently obsolete - kept for backward compability.
 942         """
 943         try:
 944             self.mlog("saved analysis %s" % self.analysisID)
 945         except:
 946             self.mlog("could not save %s" % self.analysisID)
 947
 948
 949     def deleteAnalysis(self, analysisID="default"):
 950         """ delete all of the analysis either named aName, or matching
 951             an experiment (the current one by default.) Currently a NO-OP. Obsolete.
 952         """
 953         pass
 954
 955
 956     def buildCoappear(self, motifList=[], distance="100000", strict = False):
 957         """ Builds coappear dictionary of geneIDs where motifs occur with others. Can limit to
 958             motifs in motifList (default is all motifs), and by distance (default is 100,000 base pairs.)
 959             Results are accessed through printCoappear() and saveCoappear().
 960         """
 961         occurenceList = {}
 962         self.coappear = {}
 963         processedMotifs = []
 964         motifListLen = len(motifList)
 965         if motifListLen == 0:
 966             #use all motifs
 967             motifList = self.geneToMotifKeys(thekey="mTagID")
 968             motifListLen = len(motifList)
 969
 970         for motif in motifList:
 971             if motif not in processedMotifs:
 972                 matchList = self.motifToGene(motif)
 973                 for match in matchList:
 974                     (geneID, loc) = match
 975                     if geneID not in occurenceList:
 976                         occurenceList[geneID] = []
 977
 978                     (location, sense) = loc
 979                     occurenceList[geneID].append((location, sense, motif))
 980
 981             processedMotifs.append(motif)
 982
 983         for geneID in occurenceList:
 984             occurenceList[geneID].sort()
 985             if geneID not in self.coappear:
 986                 self.coappear[geneID] = []
 987
 988             coappearing = False
 989             differentMotifs = False
 990             coappearList = []
 991             prevOccurence = occurenceList[geneID][0]
 992             del occurenceList[geneID][0]
 993             for occurence in occurenceList[geneID]:
 994                 if occurence[0] < prevOccurence[0] + distance:
 995                     coappearing = True
 996                     if occurence[2] != prevOccurence[2]:
 997                         differentMotifs = True
 998
 999                     coappearList.append(prevOccurence)
1000                 elif coappearing:
1001                     if strict:
1002                         if differentMotifs:
1003                             coappearList.append(prevOccurence)
1004                             self.coappear[geneID].append(coappearList)
1005                     else:
1006                         coappearList.append(prevOccurence)
1007                         self.coappear[geneID].append(coappearList)
1008
1009                     coappearing = False
1010                     differentMotifs = False
1011                     coappearList = []
1012
1013                 prevOccurence = occurence
1014
1015             if coappearing:
1016                 if strict:
1017                     if differentMotifs:
1018                         coappearList.append(prevOccurence)
1019                         self.coappear[geneID].append(coappearList)
1020                 else:
1021                     coappearList.append(prevOccurence)
1022                     self.coappear[geneID].append(coappearList)
1023
1024
1025     def printCoappear(self):
1026         """ prints a formatted version of the coappear dictionary built with buildCoappear()
1027         """
1028         for geneID in self.coappear:
1029             print " ===== %s =====" % str(geneID)
1030             for occurence in self.coappear[geneID]:
1031                 print str(occurence)
1032             print " ============="
1033
1034
1035     def saveCoappear(self, fileName):
1036         """ save coappear dictionary in tabular format. Returns:
1037             index, genome, locus, pos, sense, tag in tab-separated format.
1038         """
1039         index = 0
1040         outLines = []
1041         if 1:
1042             if len(fileName) < 1:
1043                 raise IOError
1044
1045             outFile = open(fileName, "a")
1046             for geneID in self.coappear:
1047                 (genome, locus) = geneID
1048                 for occurence in self.coappear[geneID]:
1049                     index += 1
1050                     coappearLine = "%d\t%s\t%s" % (index, genome, locus)
1051                     for match in occurence:
1052                         (pos, sense, tag) = match
1053                         coappearLine += "\t%s\t%s\t%s" % (pos, sense, tag)
1054
1055                     outLines.append(coappearLine + "\n")
1056
1057             outFile.writelines(outLines)
1058             outFile.close()
1059         else:
1060             self.mlog("Could not save coappear to file %s\n" % fileName)
1061
1062
1063     def sql(self, stmt, commit=""):
1064         """ executes a SQL statement and does a commit if commit is not-null. returns any results as a list.
1065         """
1066         db = sqlite.connect(self.dbName, timeout=60)
1067         sqlc = db.cursor()
1068         if self.debugSQL:
1069             print "sql: %s" % stmt
1070
1071         sqlc.execute(stmt)
1072         res = sqlc.fetchall()
1073         if commit != "":
1074             db.commit()
1075
1076         sqlc.close()
1077         db.close()
1078
1079         return res
1080
1081
1082     def batchsql(self, stmt, batch):
1083         """ executes a list of sql statements (usually inserts) stored in the list batch with a single commit.
1084         """
1085         res = []
1086         db = sqlite.connect(self.dbName, timeout=60)
1087         sqlc = db.cursor()
1088         if self.debugSQL:
1089             print "batchsql: %s" % stmt
1090             print "batchsql: %s" % str(batch)
1091
1092         sqlc.executemany(stmt, batch)
1093         db.commit()
1094         sqlc.close()
1095         db.close()
1096
1097         return res