########################################################################### # # # C O P Y R I G H T N O T I C E # # Copyright (c) 2003-10 by: # # * California Institute of Technology # # # # All Rights Reserved. # # # # Permission is hereby granted, free of charge, to any person # # obtaining a copy of this software and associated documentation files # # (the "Software"), to deal in the Software without restriction, # # including without limitation the rights to use, copy, modify, merge, # # publish, distribute, sublicense, and/or sell copies of the Software, # # and to permit persons to whom the Software is furnished to do so, # # subject to the following conditions: # # # # The above copyright notice and this permission notice shall be # # included in all copies or substantial portions of the Software. # # # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, # # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF # # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND # # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS # # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN # # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN # # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # # SOFTWARE. # ########################################################################### # # paircomp.py from cistematic.programs import Program import os, time, tempfile from cistematic.core.motif import Motif from cistematic.core import complement from os import environ if environ.get("CISTEMATIC_TEMP"): cisTemp = environ.get("CISTEMATIC_TEMP") else: cisTemp = "/tmp" tempfile.tempdir = cisTemp class Paircomp(Program): paircompPath = Program.programRoot + "/paircomp/" motifs = [] argDict = {} seq1name = "" seq2name = "" pairdir = "" def getSettings(self): return self.argDict def setSettings(self, settings): self.motifs = [] self.argDict = settings def buildInputFiles(self): """ given a Fasta file (via the constructor or the inputFile() method), grab the first two sequences and save them as two temporary sequences. """ line = " " seq1 = "" seq2 = "" inFile = open(self.inputFilePath, "r") try: while line[0] != ">": line = inFile.readline() self.seq1name = line[:-1] line = inFile.readline() while line[0] != ">": seq1 += line.strip() line = inFile.readline() self.seq2name = line[:-1] line = inFile.readline() while line[0] != ">": seq2 += line.strip() line = inFile.readline() except: if len(seq2) > 0: pass else: print "Error processing input file" return ("", "") inFile.close() self.pairdir = tempfile.mktemp() try: os.mkdir(self.pairdir) except: self.pairdir = cisTemp self.inputFile1 = self.pairdir + "/seq1.fsa" self.inputFile2 = self.pairdir + "/seq2.fsa" seq1File = open(self.inputFile1, "w") seq1File.write(self.seq1name + "\n") seq1File.write(seq1 + "\n") seq1File.close() seq2File = open(self.inputFile2, "w") seq2File.write(self.seq2name + "\n") seq2File.write(seq2 + "\n") seq2File.close() return (seq1, seq2) def buildCommand(self): self.buildInputFiles() cmd = self.paircompPath + "paircomp " + self.inputFile1 + " " + self.inputFile2 for arg in ["windowsize", "threshold"]: cmd = cmd + " " + str(self.argDict[arg]) cmd += " %s/paircomp.out" % self.pairdir print "cmd is %s" % (cmd) return cmd #set window size def setWindowSize(self, wsize): self.argDict["windowsize"] = str(wsize) # set threshold def setThreshold(self, thresh): self.argDict["threshold"] = thresh def checkDefaults(self): if "threshold" not in self.argDict: self.argDict["threshold"] = 0.9 if "windowsize" not in self.argDict: self.argDict["windowsize"] = 10 def run(self): startTime = time.time() self.checkDefaults() self.contents = os.popen(Paircomp.buildCommand(self)).readlines() outFile = open( self.pairdir + "/paircomp.out") self.contents = outFile.readlines() outFile.close() try: filenames = os.listdir(self.pairdir) for entry in filenames: os.remove(self.pairdir + "/" + entry) os.rmdir(self.pairdir) except: print "error cleaning up directory %s" % self.pairdir stopTime = time.time() print "\nThis run took %.3f - %.3f = %.3f seconds" % (startTime, stopTime, stopTime - startTime) def getWindows(self, seqNum="1"): results = [] for line in self.contents: (seq1pos, seq2pos, matches, sense) = line[:-1].split("\t") if seqNum == "1": results.append((seq1pos, seq2pos, matches, sense)) else: results.append((seq2pos, seq1pos, matches, sense)) return results def getMotifs(self): self.motifs = [] thePWM = [] motSize = self.argDict["windowsize"] matrixRow = {"A": 0, "C": 1, "G": 2, "T": 3} (seq1, seq2) = self.buildInputFiles() currentPos = -1 for line in self.contents: (seq1pos, seq2pos, matches, sense) = line[:-1].split("\t") if int(seq1pos) > currentPos: currentPos = int(seq1pos) if len(thePWM) > 0: for pos in range(len(thePWM)): thePWM[pos][matrixRow["A"]] /= seqNum thePWM[pos][matrixRow["C"]] /= seqNum thePWM[pos][matrixRow["G"]] /= seqNum thePWM[pos][matrixRow["T"]] /= seqNum self.motifs.append(Motif(self.tagID + "-PAIRCOMP-" + str(len(self.motifs) + 1), "", thePWM, [], 1.0, str(self.argDict["threshold"]))) thePWM = [] motseq1 = seq1[currentPos:currentPos+motSize] seqNum = 1.0 for pos in range(motSize): thePWM.append([0.0, 0.0, 0.0, 0.0]) NT = motseq1[pos] thePWM[pos][matrixRow[NT]] += 1 if sense == "1": motseq2 = seq2[int(seq2pos):int(seq2pos) + motSize] else: motseq2 = complement(seq2[int(seq2pos) + 1 - motSize:int(seq2pos) + 1], motSize) for pos in range(motSize): NT = motseq2[pos] thePWM[pos][matrixRow[NT]] += 1 seqNum += 1.0 if len(thePWM) > 0: for pos in range(len(thePWM)): thePWM[pos][matrixRow["A"]] /= seqNum thePWM[pos][matrixRow["C"]] /= seqNum thePWM[pos][matrixRow["G"]] /= seqNum thePWM[pos][matrixRow["T"]] /= seqNum self.motifs.append(Motif(self.tagID + "-PAIRCOMP-" + str(len(self.motifs) + 1), "", thePWM, [], 1.0, str(self.argDict["threshold"]))) return self.motifs