cistematic/experiments/experiment.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 # python parent class for experiments
  31 try:
  32     from pysqlite2 import dbapi2 as sqlite
  33 except:
  34     from sqlite3 import dbapi2 as sqlite
  35
  36 import os, time, tempfile, string
  37
  38 from cistematic.core import retrieveSeq, retrieveSeqFeatures, fasta
  39 import cistematic.core
  40 from cistematic.core.motif import matrixRow, Motif
  41 from cistematic.genomes import Genome
  42 from cistematic import version
  43 from os import environ
  44
  45 if environ.get("CISTEMATIC_TEMP"):
  46     cisTemp = environ.get("CISTEMATIC_TEMP")
  47 else:
  48     cisTemp = "/tmp"
  49
  50 tempfile.tempdir = cisTemp
  51
  52
  53 class Experiment:
  54     """ genepool = {(genomeA, geneA1): seq1, (genomeA, geneA2): seq2, (genomeB, geneB1): seq3,...}
  55         geneFeatures = ((genomeA, geneA1):[feature0, feature1, .....], (genomeB, geneB1):[feature0, feature1,...],....}
  56         programs = [(program_object1, setting1), (program_object2, setting2), ....]
  57     """
  58     genepool = {}
  59     geneFeatures = {}
  60     programs = []
  61     genepoolID = 0
  62     workdir = ""
  63     experimentID = ""
  64     experimentType = "generic"
  65     expFile = ""
  66     geneDB = ""
  67     debugSQL = False
  68
  69
  70     def __init__(self, expID, expDBFile, geneDBFile=""):
  71         self.experimentID = expID
  72         self.expFile = expDBFile
  73         self.geneDB = geneDBFile
  74         self.maskLowerCase = False
  75         self.boundToNextGene = False
  76         print "cistematic version %s" % version
  77         self.startDB()
  78
  79
  80     def __del__(self):
  81         self.savePrograms()
  82         self.saveGeneDB()
  83         self.uncacheGeneDB()
  84         self.removeWorkdir()
  85
  86
  87     def setExperimentID(self, expID):
  88         self.mlog("changing experiment ID from %s to %s" % (self.experimentID, expID))
  89         self.experimentID = expID
  90
  91
  92     def setGeneDB(self, geneDBFile=""):
  93         self.mlog("using %s as the gene database" % geneDBFile)
  94         self.geneDB = geneDBFile
  95
  96
  97     def cacheGeneDB(self, genome):
  98         try:
  99             cistematic.core.cacheGeneDB(genome)
 100             self.mlog("cached genome %s" % genome)
 101         except:
 102             self.mlog("could not cache genome %s" % genome)
 103
 104
 105     def uncacheGeneDB(self):
 106         try:
 107             cistematic.core.uncacheGeneDB()
 108         except:
 109             self.mlog("could not uncache genomes")
 110
 111
 112     def setMaskLowerCase(self, maskValue):
 113         if maskValue == True or maskValue == 1 or maskValue == "1":
 114             self.maskLowerCase = True
 115             self.setSettings("maskLowerCase", ["1"])
 116         else:
 117             self.maskLowerCase = False
 118             self.setSettings("maskLowerCase", ["0"])
 119
 120
 121     def setBoundToNextGene(self, boundValue):
 122         if boundValue == True or boundValue == 1 or boundValue == "1":
 123             self.boundToNextGene = True
 124             self.setSettings("boundToNextGene", ["1"])
 125         else:
 126             self.boundToNextGene = False
 127             self.setSettings("boundToNextGene", ["0"])
 128
 129
 130     def setSeqParameters(self, up=0, cds=1, down=0):
 131         self.upstream = up
 132         self.cds = cds
 133         self.downstream = down
 134         if cds == 0:
 135             cdsStatus = "NO "
 136         else:
 137             cdsStatus = ""
 138
 139         self.setSettings("seq_parameters",  ["%d\t%d\t%d" % (up, cds, down)])
 140         self.mlog("setting sequence retrieval parameters to %d bp upstream, %d bp downstream, and %scds" % (up, down, cdsStatus))
 141
 142
 143     def getSeqParameters(self):
 144         return (self.upstream, self.cds, self.downstream)
 145
 146
 147     def dsetLength(self):
 148         stmt="SELECT count(*) from dataset where expID = '%s' " % self.experimentID
 149         res = self.sqlexp(stmt)
 150         try:
 151             answer = int(res[0][0])
 152         except IndexError:
 153             answer = 0
 154
 155         return answer
 156
 157
 158     def resultsLength(self):
 159         stmt="SELECT count(*) from results where expID = '%s' " % self.experimentID
 160         res = self.sqlexp(stmt)
 161         try:
 162             answer = int(res[0][0])
 163         except IndexError:
 164             answer = 0
 165
 166         return answer
 167
 168
 169     def checkMotID(self, motID):
 170         answer = False
 171         stmt = "SELECT ID from motifs where expID = '%s' and mTagID = '%s' " % (self.experimentID, motID)
 172         res = self.sqlexp(stmt)
 173         if len(res) > 0:
 174             answer = True
 175
 176         return answer
 177
 178
 179     def findMotif(self,mTagID):
 180         if self.checkMotID(mTagID):
 181             return self.makeMotif(mTagID)
 182
 183         self.mlog("could not find %s" % mTagID)
 184
 185         return ""
 186
 187
 188     def makeMotif(self, motID):
 189         mPWM = []
 190         mseqs = []
 191         stmt = "SELECT motifSeq, threshold, info from motifs where expID = '%s' and mTagID = '%s' " % (self.experimentID, motID)
 192         res = self.sqlexp(stmt)
 193         entry = res[0]
 194         (seq, threshold, info) = entry
 195         stmt = "SELECT aFreq, cFreq, gFreq, tFreq from motifPWMs where expID = '%s' and mTagID = '%s' order by position" % (self.experimentID, motID)
 196         res = self.sqlexp(stmt)
 197         for entry in res:
 198             col = [0.0, 0.0, 0.0, 0.0]
 199             (aFreq, cFreq, gFreq, tFreq) = entry
 200             col[matrixRow["A"]] = aFreq
 201             col[matrixRow["C"]] = cFreq
 202             col[matrixRow["G"]] = gFreq
 203             col[matrixRow["T"]] = tFreq
 204             mPWM.append(col)
 205
 206         stmt = "SELECT sequence from motifSequences where expID = '%s' and mTagID = '%s'  and type = 'instance' " % (self.experimentID, motID)
 207         res = self.sqlexp(stmt)
 208         for entry in res:
 209             seq = entry[0]
 210             mseqs.append(seq)
 211
 212         return Motif(motID, seq, mPWM, mseqs, threshold, info)
 213
 214
 215     def saveMotif(self, mot):
 216         motID = mot.tagID
 217         motifSeq = mot.motifSeq
 218         motifInfo = mot.info
 219         motifThreshold = mot.threshold
 220         if self.checkMotID(motID):
 221             stmt = "DELETE from motifs where expID = '%s' and mTagID = '%s' " % (self.experimentID, motID)
 222             self.sqlexp(stmt, commit=True)
 223             stmt = "DELETE from motifPWMs where expID = '%s' and mTagID = '%s' " % (self.experimentID, motID)
 224             self.sqlexp(stmt, commit=True)
 225             stmt = "DELETE from motifSequences where expID = '%s' and mTagID = '%s' " % (self.experimentID, motID)
 226             self.sqlexp(stmt, commit=True)
 227
 228         stmtList = []
 229         values = "(NULL, '%s', '%s', '%s', %f, '%s')" % (self.experimentID, motID, motifSeq, motifThreshold, motifInfo)
 230         stmt = "INSERT into motifs(ID, expID, mTagID, motifSeq, threshold, info) values %s " % values
 231         self.sqlexp(stmt, commit=True)
 232         pos = 0
 233         stmt = "INSERT into motifPWMs(ID, expID, mTagID, position, aFreq, cFreq, gFreq, tFreq) values (NULL, ?, ?, ?, ?, ?, ?, ?)"
 234         for col in mot.motifPWM:
 235             aFreq = round(col[matrixRow["A"]],4)
 236             cFreq = round(col[matrixRow["C"]],4)
 237             gFreq = round(col[matrixRow["G"]],4)
 238             tFreq = abs(1.0 - aFreq - cFreq - gFreq)
 239             stmtList.append((self.experimentID, motID, pos, aFreq, cFreq, gFreq, tFreq))
 240             pos += 1
 241
 242         self.batchsqlexp(stmt, stmtList)
 243         if len(mot.sequences) > 0:
 244             stmt = "INSERT into motifSequences(ID, expID, mTagID, sequence, type, location) values (NULL, ?, ?, ?, ?, ?) "
 245             stmtList = []
 246             for seq in mot.sequences:
 247                 stmtList.append((self.experimentID, motID, seq, "instance", "-"))
 248
 249             self.batchsqlexp(stmt, stmtList)
 250
 251
 252     def exportMotifs(self, directory=".", prefix="", suffix="mot"):
 253         stmt = "SELECT distinct mTagID from results where expID = '%s'  " % (self.experimentID)
 254         res = self.sqlexp(stmt)
 255         if prefix == "":
 256             prefix = self.experimentID
 257
 258         for entry in res:
 259             motID = entry[0]
 260             mot = self.makeMotif(motID)
 261             fileName = "%s/%s-%s.%s" % (directory, prefix, motID, suffix)
 262             self.mlog("exporting %s as %s" % (motID, fileName))
 263             mot.saveMotif(fileName)
 264
 265
 266     def exportLogos(self, directory=".", prefix=""):
 267         stmt = "SELECT distinct mTagID from results where expID = '%s'  " % (self.experimentID)
 268         res = self.sqlexp(stmt)
 269         if prefix == "":
 270             prefix = self.experimentID
 271         for entry in res:
 272             motID = entry[0]
 273             mot = self.makeMotif(motID)
 274             fileName = "%s/%s-%s" % (directory, prefix, motID)
 275             self.mlog("saving logo for %s as %s.png" % (motID, fileName))
 276             mot.saveLogo(fileName)
 277
 278
 279     def appendResults(self, mot, resultsGroup="-"):
 280         motID = mot.tagID
 281         self.saveMotif(mot)
 282         stmt = "INSERT into results(ID, expID, resultsGroup, mTagID) values (NULL, '%s', '%s', '%s') " % (self.experimentID, resultsGroup, motID)
 283         self.sqlexp(stmt, "commit")
 284
 285
 286     def getGeneDB(self, geneDBFile=""):
 287         return self.geneDB
 288
 289
 290     def getLog(self):
 291         answer = []
 292         stmt = "SELECT timestamp, message from expLog where expID = '%s' order by timestamp" % self.experimentID
 293         res = self.sqlexp(stmt)
 294         for entry in res:
 295             (timestamp, message) = entry
 296             answer.append((eval(timestamp), message))
 297
 298         return answer
 299
 300
 301     def getResults(self):
 302         answer = []
 303         stmt = "SELECT distinct mTagID from results where expID = '%s'  " % (self.experimentID)
 304         res = self.sqlexp(stmt)
 305         for entry in res:
 306             motID = entry[0]
 307             answer.append(self.makeMotif(motID))
 308
 309         return answer
 310
 311
 312     def getSettings(self):
 313         answer = {}
 314         stmt = "SELECT settingName, data from settings where expID = '%s' " % (self.experimentID)
 315         res = self.sqlexp(stmt)
 316         for entry in res:
 317             (settingName, data) = entry
 318             if settingName not in answer.keys():
 319                 answer[settingName] = []
 320
 321             answer[settingName].append(data)
 322
 323         return answer
 324
 325
 326     def getSetting(self, settingName):
 327         answer = []
 328         stmt = "SELECT data from settings where expID = '%s' and settingName = '%s' " % (self.experimentID, settingName)
 329         res = self.sqlexp(stmt)
 330         for entry in res:
 331             answer.append(entry[0])
 332
 333         return answer
 334
 335
 336     def settingsHasKey(self, thekey):
 337         answer = False
 338         stmt = "SELECT distinct ID from settings where expID = '%s' and settingName = '%s' " % (self.experimentID, thekey)
 339         res = self.sqlexp(stmt)
 340         if len(res) > 0:
 341             answer = True
 342
 343         return answer
 344
 345
 346     def setSettings(self, settingName, settingList):
 347         """ insert or replace (i.e. delete previous entry) a setting with one or more setting data.
 348         """
 349         # delete existing, insert new
 350         try:
 351             stmt = "DELETE from settings where expID = '%s' and settingName = '%s'  " % (self.experimentID, settingName)
 352             res = self.sqlexp(stmt, "commit")
 353         except:
 354             pass
 355
 356         stmtList = []
 357         stmt = "INSERT into settings (ID, expID, settingName, data) values (NULL, ?, ?, ?)"
 358         for entry in settingList:
 359             stmtList.append((self.experimentID, settingName, entry))
 360
 361         res = self.batchsqlexp(stmt, stmtList)
 362
 363         return settingName
 364
 365
 366     def setSettingsID(self, settingName, data):
 367         """ return the settingsID for the inserted settingName:data pair in the settings table.
 368         """
 369         stmt = "INSERT into settings (ID, expID, settingName, data) values (NULL, '%s', '%s', \"%s\")" % (self.experimentID, settingName, data)
 370         res = self.sqlexp(stmt, "commit")
 371
 372         return res
 373
 374
 375     def getSettingsID(self, settingID):
 376         """ get a setting by settingsID  in the settings table.
 377         """
 378         answer = ""
 379         stmt = "SELECT settingName, data from  settings where expID = '%s' and ID =  %d" % (self.experimentID, int(settingID))
 380         res = self.sqlexp(stmt)
 381         try:
 382             (name, data) = res[0]
 383             answer = (name, data)
 384         except:
 385             pass
 386
 387         return answer
 388
 389
 390     def setRun(self, progName, datasetID, settingsID):
 391         values = "(NULL, '%s', '%s', '%s', %d, '%s', '%s')" % (self.experimentID, progName, datasetID, settingsID, time.localtime(), "-")
 392         stmt = "INSERT into runs (ID, expID, progName, datasetGroup, settingsID, timestamp, resultsGroup) values %s" % values
 393         runID = self.sqlexp(stmt, "commit")
 394         self.mlog("run %s: program %s with settings %d and dataset %s" % (runID, progName, settingsID, datasetID))
 395
 396         return int(runID)
 397
 398
 399     def getRun(self, rID):
 400         res = ""
 401         stmt = "SELECT progName, datasetGroup, settingsID, timestamp, resultsGroup from runs where expID = '%s' and ID = %d" % (self.experimentID, rID)
 402         res = self.sqlexp(stmt)
 403         (progName, datasetID, settingsID, timestamp, resultsID) = res[0]
 404         if resultsID == "-":
 405             resultsID = []
 406             if datasetID != "chromolist":
 407                 datasetID = int(datasetID)
 408
 409         return (progName, datasetID, int(settingsID), timestamp, resultsID)
 410
 411
 412     def getRunsByProg(self, prog):
 413         runs = self.getRuns()
 414         matchingRuns = []
 415         for entry in runs:
 416             if (runs[entry][0] == prog):
 417                 matchingRuns.append(runs[entry])
 418
 419         return matchingRuns
 420
 421
 422     def getRuns(self):
 423         runs = {}
 424         res = ''
 425         stmt = "SELECT ID, progName, datasetGroup, settingsID, timestamp, resultsGroup from runs where expID = '%s' " % (self.experimentID)
 426         res = self.sqlexp(stmt)
 427         for entry in res:
 428             (ID, progName, datasetID, settingsID, timestamp, resultsID) = entry
 429             runs[int(ID)] = (progName, datasetID, settingsID, eval(timestamp), resultsID)
 430
 431         return runs
 432
 433
 434     def appendDataset(self, dset):
 435         stmtList = []
 436         stmt = "INSERT into dataset(ID, genome, locus, sequence, expID) values (NULL, ?, ?, ?, ?)  "
 437         for (genome, entries) in dset:
 438             for entry in entries:
 439                 if len(entry) == 2:
 440                     (locus, sequence) = entry
 441                 else:
 442                     locus = entry
 443                     sequence = "-"
 444
 445                 stmtList.append((genome, locus, sequence, self.experimentID))
 446
 447         res = self.batchsqlexp(stmt, stmtList)
 448
 449
 450     def getDataset(self):
 451         dset = []
 452         stmt = "SELECT genome, locus, sequence from dataset where expID = '%s' " % (self.experimentID)
 453         res = self.sqlexp(stmt)
 454         for entry in res:
 455             (genome, locus, sequence) = entry
 456             if sequence == "-":
 457                 dset.append((genome, [locus]))
 458             else:
 459                 dset.append((genome, [locus, sequence]))
 460
 461         return dset
 462
 463
 464     def getDatasetNames(self):
 465         dset = []
 466         stmt = "select settingName from settings where settingName like 'dataset" + "%' and expID = '" + self.experimentID + "' "
 467         res = self.sqlexp(stmt)
 468         for entry in res:
 469             dset.append(entry)
 470
 471         return dset
 472
 473
 474     def getDatasetIDs(self):
 475         dset = []
 476         stmt = "select ID from settings where settingName like 'dataset" + "%' and expID = '" + self.experimentID + "' "
 477         res = self.sqlexp(stmt)
 478         for entry in res:
 479             dset.append(entry)
 480
 481         return dset
 482
 483
 484     def getFeatures(self, geneID):
 485         results = []
 486         try:
 487             res = []
 488             stmt = "select featureType, start, stop, orientation from sequenceFeatures where seqGenome = '%s' and seqID = '%s' and expID = '%s' " % (geneID[0], geneID[1], self.experimentID)
 489             res = self.sqlexp(stmt)
 490             for entry in res:
 491                 (ftype, start, stop, orientation) = entry
 492                 results.append((ftype, start, stop, orientation))
 493         except:
 494             pass
 495
 496         results.sort()
 497
 498         return results
 499
 500
 501     def absoluteLocation(self, match, featureLength, seqparams="", gidCoordinates="", customDB=""):
 502         """ Returns the absolute location of the start of a match = ((genome, gID), (pos, sense)), given a relative
 503             location with respect to the gene and a feature length.
 504             Can be passed cached sequence parameters (up, cds, down) and geneEntry tuple to avoid hitting database.
 505         """
 506         result = ["", 0, "F"]
 507         (geneID, loc) = match
 508         (pos, sense) = loc
 509         pos = int(pos)
 510         if seqparams == "":
 511             seqparams = self.getSeqParameters()
 512
 513         (up, cds, down) = seqparams
 514         if gidCoordinates == "":
 515             gidCoordinates = cistematic.core.geneEntry(geneID)
 516
 517         (gidChrom, gidStart, gidStop, gidLength, gidSense) = gidCoordinates
 518         result[0] = gidChrom
 519         if self.boundToNextGene:
 520             up = cistematic.core.upstreamToNextGene(geneID, up, db=customDB)
 521             down = cistematic.core.downstreamToNextGene(geneID, down, db=customDB)
 522
 523         if gidSense == "F":
 524             result[2] = sense
 525             if pos < up or cds > 0:
 526                 result[1] = gidStart - up + pos
 527             else:
 528                 result[1] = gidStop - up + pos
 529         else:
 530             if gidStart > gidStop:
 531                 temp = gidStart
 532                 gidStart = gidStop
 533                 gidStop = temp
 534
 535             if sense == "F":
 536                 result[2] = "R"
 537
 538             if pos < up or cds > 0:
 539                 result[1] = gidStop + up - pos - featureLength
 540             else:
 541                 result[1] = gidStart + up - pos - featureLength
 542
 543         return result
 544
 545
 546     def setWorkdir(self, wdir=""):
 547         self.workdir = wdir
 548         self.mlog("changed workdir to %s" % (wdir))
 549
 550
 551     def resetDataset(self):
 552         stmt = "DELETE from dataset where expID = '%s'" % self.experimentID
 553         res = self.sqlexp(stmt, "commit")
 554         self.mlog("log reset")
 555
 556
 557     def resetLog(self):
 558         stmt = "DELETE from expLog where expID = '%s'" % self.experimentID
 559         res = self.sqlexp(stmt, "commit")
 560         self.mlog("log reset")
 561
 562
 563     def resetResults(self):
 564         stmt = "DELETE from results where expID = '%s'" % self.experimentID
 565         res = self.sqlexp(stmt, "commit")
 566         self.mlog("results reset")
 567
 568
 569     def resetRuns(self):
 570         stmt = "DELETE from expLog where expID = '%s'" % self.experimentID
 571         res = self.sqlexp(stmt, "commit")
 572         self.mlog("runs reset")
 573
 574
 575     def resetPrograms(self):
 576         self.programs = []
 577         self.mlog("programs reset")
 578
 579
 580     def resetSettings(self):
 581         stmt = "DELETE from expLog where expID = '%s'" % self.experimentID
 582         res = self.sqlexp(stmt, "commit")
 583         self.log("settings reset")
 584
 585
 586     def appendProgram(self, program):
 587         self.programs.append((program, self.setSettingsID(program.name(), program.getSettings())))
 588         self.mlog("adding program %s" % program.name())
 589
 590
 591     def removeProgram(self, progName):
 592         for index in range(len(self.programs)):
 593             if self.programs[index][0].name() == progName:
 594                 del self.programs[index]
 595                 self.mlog("removed program %s" % progName)
 596
 597
 598     def createWorkdir(self):
 599         if self.workdir == "":
 600             self.workdir = tempfile.mktemp()
 601
 602         try:
 603             os.mkdir(self.workdir)
 604             self.mlog("created workdir: %s" % (self.workdir))
 605         except:
 606             pass
 607
 608
 609     def removeWorkdir(self):
 610         try:
 611             filenames = os.listdir(self.workdir)
 612             for entry in filenames:
 613                 os.remove("%s/%s" % (self.workdir, entry))
 614             os.rmdir(self.workdir)
 615         except:
 616             self.mlog("could not delete workdir: %s" % (self.workdir))
 617
 618
 619     def mlog(self,msg):
 620         print msg
 621         stmt = "INSERT into expLog(ID, expID, timestamp, message) values (NULL, '%s', '%s', '%s') " % (self.experimentID, time.localtime(),string.replace(msg, "'", '"'))
 622         self.sqlexp(stmt, "commit")
 623
 624
 625     def logToString(self):
 626         response = ""
 627         for line in self.getLog():
 628             response +=  time.asctime(line[0]) + ": " + line[1] + "\n"
 629
 630         return response
 631
 632
 633     def printLog(self):
 634         for line in self.getLog():
 635             print "%s: %s" % (time.asctime(line[0]), line[1])
 636
 637
 638     def tailLog(self):
 639         theLog = self.getLog()
 640         for line in theLog[-10:]:
 641             print "%s: %s" % (time.asctime(line[0]), line[1])
 642
 643
 644     def loadPrograms(self):
 645         self.programs = []
 646         for progs in self.getSetting("loaded_programs"):
 647             (progs0, progs1, progs2) = progs.split("\t")
 648             execString = "from %s import %s" % (progs0, progs1)
 649             exec execString
 650             execString = 'self.programs.append((apply(%s), %s))' % (progs1, progs2)
 651             exec execString
 652
 653
 654     def savePrograms(self):
 655         progs = []
 656         for (program, settingID) in self.programs:
 657             progs.append("%s\t%s\t%s" % (program.__class__.__module__, program.__class__.__name__, settingID))
 658
 659         if len(progs) > 0:
 660             self.setSettings("loaded_programs", progs)
 661
 662     def saveGeneDB(self):
 663         if self.geneDB != "":
 664             self.setSettings("geneDB", [self.geneDB])
 665
 666
 667     def loadFasta(self, fastaFile, genomeName):
 668         """ load fasta file into the dataset.
 669         """
 670         genIDList = self.loadFastaFromFile(fastaFile, genomeName)
 671
 672         return [(genomeName, genIDList)]
 673
 674
 675     def loadFastaFromFile(self, fastaFile, genomeName):
 676         seqArray = []
 677         seqLen = 0
 678         seqName = []
 679         geneDBPath = "%s.genedb" % genomeName
 680         self.mlog("Loading fasta file %s into database %s" % (genomeName, geneDBPath))
 681         aGenome = Genome(genomeName, dbFile=geneDBPath)
 682         aGenome.createGeneDB()
 683         inFile = open(fastaFile, "r")
 684         header = inFile.readline()
 685         while header != "":
 686             seqArray = []
 687             seqLen = 0
 688             chromID = header.strip()[1:]
 689             currentLine = inFile.readline()
 690             while currentLine != "" and currentLine[0] != ">":
 691                 lineSeq = currentLine.strip()
 692                 seqLen += len(lineSeq)
 693                 seqArray.append(lineSeq)
 694                 currentLine = inFile.readline()
 695
 696             seq = string.join(seqArray, "")
 697             print "Added sequence %s to database" % chromID
 698             aGenome.addSequence((genomeName, chromID), seq, "chromosome", str(seqLen))
 699             aGenome.addChromosomeEntry(chromID, chromID, "db")
 700             aGenome.addGeneEntry((genomeName, chromID), chromID, 0, seqLen - 1, "F", "IMPORT", "1")
 701             header = currentLine
 702             seqName.append(chromID)
 703
 704         inFile.close()
 705         aGenome.createIndices()
 706         self.setGeneDB(geneDBPath)
 707
 708         return seqName
 709
 710
 711     def loadGenepool(self):
 712         self.genepool = {}
 713         for geneTuple in self.getDataset():
 714             genome = geneTuple[0]
 715             geneList = geneTuple[1]
 716             for gene in geneList:
 717                 try:
 718                     if len(gene) == 2:
 719                         tag = gene[0]
 720                         geneID = (genome, tag)
 721                         seq = gene[1]
 722                         # need to deal with masking lowercase
 723                         self.genepool[geneID] = seq.upper()
 724                     else:
 725                         geneID = (genome, gene)
 726                         # note that we will only keep one copy of a geneid in the genepool, even if we
 727                         # retrieve it multiple times.
 728                         self.genepool[geneID] = retrieveSeq(geneID, self.upstream, self.cds, self.downstream, self.geneDB, self.maskLowerCase, self.boundToNextGene)
 729                 except:
 730                     self.mlog("could not load %s" % (str(geneID)))
 731
 732         self.genepoolID = self.setSettingsID("genepool", self.genepool.keys())
 733
 734
 735     def loadGeneFeatures(self):
 736         stmt = "DELETE from sequenceFeatures where expID = '%s'" % self.experimentID
 737         res = self.sqlexp(stmt, "commit")
 738         stmtList = []
 739         stmt = "INSERT into sequenceFeatures (ID, expID, seqGenome, seqID, featureType, start, stop, orientation) values (NULL, ?, ?, ?, ?, ?, ?, ?)"
 740         for geneTuple in self.getDataset():
 741             genome = geneTuple[0]
 742             geneList = geneTuple[1]
 743             for gene in geneList:
 744                 try:
 745                     geneID = (genome, gene)
 746                     results = retrieveSeqFeatures(geneID, self.upstream, self.cds, self.downstream, self.boundToNextGene, self.geneDB)
 747                     for entry in results:
 748                         (ftype, start, stop, orientation) = entry
 749                         stmtList.append((self.experimentID, genome, gene, ftype, start, stop, orientation))
 750                 except:
 751                     self.mlog("could not find features for %s" % (gene))
 752
 753         if len(stmtList) > 0:
 754             self.batchsqlexp(stmt, stmtList)
 755
 756
 757     def reload(self):
 758         pass
 759
 760
 761     def toFile(self, geneIDList, filename, geneDict=[]):
 762         outFile = open(filename, "w")
 763         for geneID in geneIDList:
 764             if geneID in geneDict:
 765                 outFile.write(fasta(geneID, geneDict[geneID]))
 766             elif geneID in self.genepool:
 767                 outFile.write(fasta(geneID, self.genepool[geneID]))
 768             else:
 769                 self.mlog("could not write %s to file" % str(geneID))
 770
 771         outFile.close()
 772
 773
 774     def createDataFile(self, datasetID=-1, geneIDList=[], geneDict=[]):
 775         oldtempdir = tempfile.tempdir
 776         tempfile.tempdir = self.workdir
 777         dataFile = tempfile.mktemp()
 778         tempfile.tempdir = oldtempdir
 779         if datasetID < 0:
 780             datasetID = self.genepoolID
 781
 782         if len(geneIDList) < 1:
 783             settingsList =  self.getSettingsID(datasetID)
 784             geneIDList = eval(settingsList[1])
 785
 786         try:
 787             self.toFile(geneIDList, dataFile, geneDict)
 788         except:
 789             self.mlog("could not create dataFile %s" % (dataFile))
 790
 791         return dataFile
 792
 793
 794     def initialize(self, dataset=[], workdir=""):
 795         if workdir != "":
 796             self.setWorkdir(workdir)
 797
 798         self.createWorkdir()
 799         if len(dataset) > 0:
 800             self.appendDataset(dataset)
 801             self.loadGenepool()
 802             self.loadGeneFeatures()
 803
 804
 805     def run(self):
 806         if len(self.programs)  == 0:
 807             self.mlog("Must instantiate one or more programs first")
 808         elif len(self.genepool) == 0:
 809             self.mlog("Must have one or more valid sequences in the dataset")
 810
 811
 812     def sqlexp(self, stmt, commit=""):
 813         db = sqlite.connect(self.expFile, timeout=60)
 814         sqlc = db.cursor()
 815         if self.debugSQL:
 816             print "sqlexp: %s" % stmt
 817
 818         sqlc.execute(stmt)
 819         res = sqlc.fetchall()
 820         if commit != "":
 821             db.commit()
 822
 823         if stmt[0:6] == "INSERT":
 824             res = sqlc.lastrowid
 825
 826         sqlc.close()
 827         db.close()
 828
 829         return res
 830
 831
 832     def batchsqlexp(self, stmt, batch):
 833         """ executes a list of sql statements (usually inserts) stored in the list batch with a single commit.
 834         """
 835         res = []
 836         db = sqlite.connect(self.expFile, timeout=60)
 837         sqlc = db.cursor()
 838         if self.debugSQL:
 839             print "batchsql: %s" % stmt
 840             print "batchsql: %s" % str(batch)
 841
 842         sqlc.executemany(stmt, batch)
 843         db.commit()
 844         sqlc.close()
 845         db.close()
 846
 847         return res
 848
 849
 850     def setExternalStatus(self, status):
 851         self.mlog("Setting status to %s" % status)
 852         statfile = open("%s.status" % self.expFile, "w")
 853         statfile.write(status)
 854         statfile.close()
 855
 856
 857     def resetExternalStatus(self):
 858         try:
 859             os.remove("%s.status" % self.expFile)
 860         except:
 861             pass
 862
 863
 864     def startDB(self):
 865         if not os.path.exists(self.expFile):
 866             try:
 867                 db = sqlite.connect(self.expFile, timeout=60)
 868                 sql = db.cursor()
 869                 sql.execute("CREATE table experiment(ID INTEGER PRIMARY KEY, expID varchar, expType varchar, expStatus varchar, timestamp varchar)")
 870                 sql.execute("CREATE table dataset(ID INTEGER PRIMARY KEY, expID varchar, datasetGroup varchar, genome varchar, locus varchar, sequence varchar)")
 871                 sql.execute("CREATE table results(ID INTEGER PRIMARY KEY, expID varchar, resultsGroup varchar, mTagID varchar)")
 872                 sql.execute("CREATE table motifs(ID INTEGER PRIMARY KEY, expID varchar, mTagID varchar, motifSeq varchar, threshold varchar, info varchar)")
 873                 sql.execute("CREATE table motifPWMs(ID INTEGER PRIMARY KEY, expID varchar, mTagID varchar, position int, aFreq float, cFreq float, gFreq float, tFreq float)")
 874                 sql.execute("CREATE table motifSequences(ID INTEGER PRIMARY KEY, expID varchar, mTagID varchar, sequence varchar, type varchar, location varchar)")
 875                 sql.execute("CREATE table settings(ID INTEGER PRIMARY KEY, expID varchar, settingName varchar, data varchar)")
 876                 sql.execute("CREATE table runs(ID INTEGER PRIMARY KEY, expID varchar, progName varchar, datasetGroup varchar, settingsID int, timestamp varchar, resultsGroup varchar)")
 877                 sql.execute("CREATE table expLog(ID INTEGER PRIMARY KEY, expID varchar, timestamp varchar, message varchar)")
 878                 sql.execute("CREATE table sequenceFeatures(ID INTEGER PRIMARY KEY, expID varchar, seqGenome varchar, seqID varchar, featureType varchar, start int, stop int, orientation varchar)")
 879
 880                 sql.execute("CREATE INDEX settingIndex1 on settings(expID, settingName)")
 881                 sql.execute("CREATE INDEX datasetIndex1 on dataset(expID, datasetGroup)")
 882                 sql.execute("CREATE INDEX motifsIndex1 on motifs(expID, mTagID)")
 883                 sql.execute("CREATE INDEX motifPWMsIndex1 on motifPWMs(expID, mTagID)")
 884                 sql.execute("CREATE INDEX motifSequencesIndex1 on motifSequences(expID, mTagID)")
 885                 sql.execute("CREATE INDEX featuresIndex1 on sequenceFeatures(expID, seqGenome, seqID)")
 886
 887                 sql.execute("INSERT INTO settings(ID, expID, settingName, data) values (NULL, '%s', 'experimentType', '%s')" % (self.experimentID, self.experimentType))
 888
 889                 db.commit()
 890                 sql.close()
 891                 db.close()
 892                 self.mlog("Created experiment database %s" % self.expFile)
 893             except:
 894                 self.mlog("Could not create experiment database %s" % self.expFile)
 895         else:
 896             self.mlog("Using existing experiment database %s" % self.expFile)
 897
 898         if self.settingsHasKey("loaded_programs"):
 899             self.loadPrograms()
 900
 901         if self.settingsHasKey("seq_parameters"):
 902             res = self.getSetting("seq_parameters")
 903             (up, cds, down) = res[0].split("\t")
 904             self.setSeqParameters(int(up), int(cds), int(down))
 905         else:
 906             self.setSeqParameters()
 907
 908         if self.settingsHasKey("maskLowerCase"):
 909             res = self.getSetting("maskLowerCase")
 910             self.setMaskLowerCase(res[0])
 911         else:
 912             self.setMaskLowerCase(False)
 913
 914         if self.settingsHasKey("boundToNextGene"):
 915             res = self.getSetting("boundToNextGene")
 916             self.setBoundToNextGene(res[0])
 917         else:
 918             self.setBoundToNextGene(False)
 919
 920         if self.settingsHasKey("experimentType"):
 921             res = self.getSetting("experimentType")
 922             self.experimentType = res[0]
 923
 924         if self.settingsHasKey("geneDB"):
 925             res = self.getSetting("geneDB")
 926             self.geneDB = res[0]
 927
 928         if self.dsetLength() > 0:
 929             self.loadGenepool()
 930
 931         self.createWorkdir()