Mercurial > repos > yufei-luo > s_mart

comparison SMART/Java/Python/GetFlanking.py @ 46:169d364ddd91

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Mon, 30 Sep 2013 03:19:26 -0400
parents 44d5973c188c
children
comparison
equal deleted inserted replaced
45:e454402ba9d9 46:169d364ddd91
42 TAG_DISTANCE = "distance_" 42 TAG_DISTANCE = "distance_"
43 TAG_SENSE = "_sense" 43 TAG_SENSE = "_sense"
44 TAG_REGION = "_region" 44 TAG_REGION = "_region"
45 TAGS_REGION = {-1: "_upstream", 0: "", 1: "_downstream"} 45 TAGS_REGION = {-1: "_upstream", 0: "", 1: "_downstream"}
46 TAGS_RREGION = {-1: "upstream", 0: "overlapping", 1: "downstream"} 46 TAGS_RREGION = {-1: "upstream", 0: "overlapping", 1: "downstream"}
47 TAGS_SENSE = {-1: "antisense", 0: "", 1: "collinear"} 47 TAGS_SENSE = {-1: "antisense", 0: "", 1: "colinear"}
48 STRANDSTOSTR = {-1: "(-)", 0: "", 1: "(+)"} 48 STRANDSTOSTR = {-1: "(-)", 0: "", 1: "(+)"}
49 49
50 50
51 def getOrderKey(transcript, direction, input): 51 def getOrderKey(transcript, direction):
52 if direction == 1: 52 if direction == 1:
53 if input == QUERY: 53 return transcript.getEnd()
54 return (transcript.getEnd(), -transcript.getStart()) 54 return - transcript.getStart()
55 return (transcript.getStart(), -transcript.getEnd()) 55
56 if input == QUERY: 56 def isInGoodRegion(transcriptRef, transcriptQuery, direction):
57 return (-transcript.getStart(), transcript.getEnd()) 57 if direction == 1:
58 return (-transcript.getEnd(), transcript.getStart()) 58 return transcriptQuery.getEnd() > transcriptRef.getEnd()
59 return transcriptQuery.getStart() < transcriptRef.getStart()
59 60
60 61
61 class GetFlanking(object): 62 class GetFlanking(object):
62 63
63 def __init__(self, verbosity): 64 def __init__(self, verbosity):
64 self.verbosity = verbosity 65 self.verbosity = verbosity
65 self.transcripts = dict([id, {}] for id in INPUTS) 66 self.transcripts = dict([id, {}] for id in INPUTS)
66 self.directions = [] 67 self.directions = []
67 self.noOverlap = False 68 self.noOverlap = False
68 self.colinear = False 69 self.colinear = False
69 self.antisense = False 70 self.antisense = False
70 self.distance = None 71 self.distance = None
71 self.minDistance = None 72 self.minDistance = None
72 self.maxDistance = None 73 self.maxDistance = None
73 self.tagName = "flanking" 74 self.tagName = "flanking"
74 75
75 def setInputFile(self, fileName, format, id): 76 def setInputFile(self, fileName, format, id):
76 chooser = ParserChooser(self.verbosity) 77 chooser = ParserChooser(self.verbosity)
77 chooser.findFormat(format) 78 chooser.findFormat(format)
78 parser = chooser.getParser(fileName) 79 parser = chooser.getParser(fileName)
79 for transcript in parser.getIterator(): 80 for transcript in parser.getIterator():
80 chromosome = transcript.getChromosome() 81 chromosome = transcript.getChromosome()
81 if chromosome not in self.transcripts[id]: 82 if chromosome not in self.transcripts[id]:
82 self.transcripts[id][chromosome] = [] 83 self.transcripts[id][chromosome] = []
83 self.transcripts[id][chromosome].append(transcript) 84 self.transcripts[id][chromosome].append(transcript)
84 85
85 def setOutputFile(self, fileName): 86 def setOutputFile(self, fileName):
86 self.writer = TranscriptWriter(fileName, "gff3", self.verbosity) 87 self.writer = TranscriptWriter(fileName, "gff3", self.verbosity)
87 88
88 def addUpstreamDirection(self, upstream): 89 def addUpstreamDirection(self, upstream):
89 if upstream: 90 if upstream:
90 self.directions.append(-1) 91 self.directions.append(-1)
91 92
92 def addDownstreamDirection(self, downstream): 93 def addDownstreamDirection(self, downstream):
93 if downstream: 94 if downstream:
94 self.directions.append(1) 95 self.directions.append(1)
95 96
96 def setColinear(self, colinear): 97 def setColinear(self, colinear):
97 self.colinear = colinear 98 self.colinear = colinear
98 99
99 def setAntisense(self, antisense): 100 def setAntisense(self, antisense):
100 self.antisense = antisense 101 self.antisense = antisense
101 102
102 def setNoOverlap(self, noOverlap): 103 def setNoOverlap(self, noOverlap):
103 self.noOverlap = noOverlap 104 self.noOverlap = noOverlap
104 105
105 def setMinDistance(self, distance): 106 def setMinDistance(self, distance):
106 self.minDistance = distance 107 self.minDistance = distance
107 108
108 def setMaxDistance(self, distance): 109 def setMaxDistance(self, distance):
109 self.maxDistance = distance 110 self.maxDistance = distance
110 111
111 def setNewTagName(self, tagName): 112 def setNewTagName(self, tagName):
112 self.tagName = tagName 113 self.tagName = tagName
113 114
114 def match(self, transcriptQuery, transcriptRef, direction): 115 def match(self, transcriptRef, transcriptQuery, direction):
115 #print "comparing", transcriptQuery, "with", transcriptRef, "on direction", direction 116 if self.noOverlap and transcriptRef.overlapWith(transcriptQuery):
116 if direction == 1 and transcriptRef.getEnd() < transcriptQuery.getStart(): 117 return False
117 return False 118 if self.colinear and transcriptRef.getDirection() != transcriptQuery.getDirection():
118 if direction == -1 and transcriptQuery.getEnd() < transcriptRef.getStart(): 119 return False
119 return False 120 if self.antisense and transcriptRef.getDirection() == transcriptQuery.getDirection():
120 if self.noOverlap and transcriptRef.overlapWith(transcriptQuery): 121 return False
121 return False 122 if self.minDistance != None or self.maxDistance != None:
122 if self.colinear and transcriptRef.getDirection() != transcriptQuery.getDirection(): 123 distance = transcriptRef.getDistance(transcriptQuery)
123 return False 124 if self.minDistance != None and distance < self.minDistance:
124 if self.antisense and transcriptRef.getDirection() == transcriptQuery.getDirection(): 125 return False
125 return False 126 if self.maxDistance != None and distance > self.maxDistance:
126 if self.minDistance != None or self.maxDistance != None: 127 return False
127 distance = transcriptRef.getDistance(transcriptQuery) 128 return True
128 if self.minDistance != None and distance < self.minDistance: 129
129 return False 130 def getFlanking(self, direction):
130 if self.maxDistance != None and distance > self.maxDistance: 131 for chromosome in sorted(self.transcripts[REFERENCE].keys()):
131 return False 132 if chromosome not in self.transcripts[QUERY]:
132 return True 133 continue
133 134 sortedTranscripts = dict([id, {}] for id in INPUTS)
134 def getFlanking(self, chromosome, direction): 135 for id in INPUTS:
135 if chromosome not in self.transcripts[REFERENCE]: 136 sortedTranscripts[id] = sorted(self.transcripts[id][chromosome], key = lambda t: getOrderKey(t, direction))
136 return 137 refIndex = 0
137 sortedTranscripts = dict([id, {}] for id in INPUTS) 138 currentRefs = []
138 for id in INPUTS: 139 outputs = set()
139 sortedTranscripts[id] = sorted(self.transcripts[id][chromosome], key = lambda t: getOrderKey(t, direction, id)) 140 progress = Progress(len(sortedTranscripts[QUERY]), "Reading chr %s %s" % (chromosome, STRANDSTOSTR[direction]), self.verbosity)
140 refIndex = 0 141 for query in sortedTranscripts[QUERY]:
141 progress = Progress(len(sortedTranscripts[QUERY]), "Reading chr %s %s" % (chromosome, STRANDSTOSTR[direction]), self.verbosity) 142 while refIndex < len(sortedTranscripts[REFERENCE]) and isInGoodRegion(sortedTranscripts[REFERENCE][refIndex], query, direction):
142 for query in sortedTranscripts[QUERY]: 143 currentRefs.append(sortedTranscripts[REFERENCE][refIndex])
143 #print "Q: ", query 144 refIndex += 1
144 #print "R1: ", sortedTranscripts[REFERENCE][refIndex] 145 nextCurrentRefs = []
145 while not self.match(query, sortedTranscripts[REFERENCE][refIndex], direction): 146 for currentRef in currentRefs:
146 refIndex += 1 147 if self.match(currentRef, query, direction):
147 if refIndex == len(sortedTranscripts[REFERENCE]): 148 if currentRef not in self.flankings:
148 progress.done() 149 self.flankings[currentRef] = {}
149 #print "done" 150 self.flankings[currentRef][direction * currentRef.getDirection()] = query
150 return 151 else:
151 #print "R2: ", sortedTranscripts[REFERENCE][refIndex] 152 nextCurrentRefs.append(currentRef)
152 self.flankings[query][direction] = sortedTranscripts[REFERENCE][refIndex] 153 currentRefs = nextCurrentRefs
153 progress.inc() 154 progress.inc()
154 progress.done() 155 progress.done()
155 156
156 def setTags(self, query, reference, direction): 157 def setTags(self, query, reference, direction):
157 refName = reference.getTagValue("ID") 158 refName = reference.getTagValue("ID")
158 if refName == None: 159 if refName == None:
159 refName = reference.getName() 160 refName = reference.getName()
160 if refName == None: 161 if refName == None:
161 refName = reference.__str__() 162 refName = reference.__str__()
162 query.setTagValue("%s%s" % (self.tagName, TAGS_REGION[direction*query.getDirection()]), refName) 163 query.setTagValue("%s%s" % (self.tagName, TAGS_REGION[direction]), refName)
163 query.setTagValue("%s_%s%s" % (TAG_DISTANCE, self.tagName, TAGS_REGION[direction*query.getDirection()]), query.getDistance(reference)) 164 query.setTagValue("%s_%s%s" % (TAG_DISTANCE, self.tagName, TAGS_REGION[direction]), query.getDistance(reference))
164 query.setTagValue("%s_%s" % (TAG_SENSE, self.tagName), TAGS_SENSE[query.getDirection() * reference.getDirection()]) 165 if direction == 0:
165 if direction == 0: 166 query.setTagValue("%s_%s" % (TAG_SENSE, self.tagName), TAGS_SENSE[query.getDirection() * reference.getDirection()])
166 query.setTagValue("%s_%s" % (TAG_REGION, self.tagName), TAGS_RREGION[cmp(query.getRelativeDistance(reference), 0)]) 167 query.setTagValue("%s_%s" % (TAG_REGION, self.tagName), TAGS_RREGION[cmp(query.getRelativeDistance(reference), 0)])
167 for tag in reference.getTagNames(): 168 for tag in reference.getTagNames():
168 if tag not in ("quality", "feature"): 169 if tag not in ("quality", "feature"):
169 query.setTagValue("%s%s_%s" % (self.tagName, TAGS_REGION[direction*query.getDirection()], tag), reference.getTagValue(tag)) 170 query.setTagValue("%s%s_%s" % (self.tagName, TAGS_REGION[direction], tag), reference.getTagValue(tag))
170 return query 171 return query
171 172
172 def write(self): 173 def write(self):
173 progress = Progress(len(self.flankings.keys()), "Printing data", self.verbosity) 174 outputs = set()
174 for transcriptQuery in self.flankings.keys(): 175 progress = Progress(len(self.flankings.keys()), "Printing data", self.verbosity)
175 if not self.flankings[transcriptQuery]: 176 for transcriptRef in self.flankings.keys():
176 self.writer.addTranscript(transcriptQuery) 177 if self.directions:
177 elif self.directions: 178 for direction in self.directions:
178 for direction in self.directions: 179 if direction in self.flankings[transcriptRef]:
179 #relativeDirection = direction if transcriptQuery.getDirection() == 1 else - direction 180 query = self.flankings[transcriptRef][direction]
180 relativeDirection = direction * transcriptQuery.getDirection() 181 outputs.add(self.setTags(query, transcriptRef, direction))
181 if relativeDirection in self.flankings[transcriptQuery]: 182 else:
182 transcriptRef = self.flankings[transcriptQuery][relativeDirection] 183 if self.flankings[transcriptRef]:
183 transcriptQuery = self.setTags(transcriptQuery, transcriptRef, relativeDirection) 184 query = sorted(self.flankings[transcriptRef].values(), key = lambda query: query.getDistance(transcriptRef))[0]
184 self.writer.addTranscript(transcriptQuery) 185 outputs.add(self.setTags(query, transcriptRef, 0))
185 else: 186 progress.inc()
186 transcriptRef = sorted(self.flankings[transcriptQuery].values(), key = lambda transcriptRef: transcriptQuery.getDistance(transcriptRef))[0] 187 for transcript in sorted(list(outputs), key = lambda flanking: (flanking.getChromosome(), flanking.getStart(), flanking.getEnd())):
187 self.writer.addTranscript(self.setTags(transcriptQuery, transcriptRef, 0)) 188 self.writer.addTranscript(transcript)
188 progress.inc() 189 self.writer.close()
189 progress.done() 190 progress.done()
190 191
191 def run(self): 192 def run(self):
192 for chromosome in sorted(self.transcripts[QUERY].keys()): 193 self.flankings = {}
193 self.flankings = dict([query, {}] for query in self.transcripts[QUERY][chromosome]) 194 for direction in STRANDS:
194 for direction in STRANDS: 195 self.getFlanking(direction)
195 #print "comparison", chromosome, direction 196 self.write()
196 self.getFlanking(chromosome, direction)
197 self.write()
198 self.writer.close()
199 197
200 if __name__ == "__main__": 198 if __name__ == "__main__":
201 199
202 description = "Get Flanking v1.0.1: Get the flanking regions of a set of reference. [Category: Data Selection]" 200 description = "Get Flanking v1.0.1: Get the flanking regions of a set of reference. [Category: Data Selection]"
203 201
204 parser = OptionParser(description = description) 202 parser = OptionParser(description = description)
205 parser.add_option("-i", "--input1", dest="inputFileName1", action="store", type="string", help="query input file [compulsory] [format: file in transcript format given by -f]") 203 parser.add_option("-i", "--input1", dest="inputFileName1", action="store", type="string", help="query input file [compulsory] [format: file in transcript format given by -f]")
206 parser.add_option("-f", "--format1", dest="format1", action="store", type="string", help="format of previous file [compulsory] [format: transcript file format]") 204 parser.add_option("-f", "--format1", dest="format1", action="store", type="string", help="format of previous file [compulsory] [format: transcript file format]")
207 parser.add_option("-j", "--input2", dest="inputFileName2", action="store", type="string", help="reference input file [compulsory] [format: file in transcript format given by -g]") 205 parser.add_option("-j", "--input2", dest="inputFileName2", action="store", type="string", help="reference input file [compulsory] [format: file in transcript format given by -g]")
208 parser.add_option("-g", "--format2", dest="format2", action="store", type="string", help="format of previous file [compulsory] [format: transcript file format]") 206 parser.add_option("-g", "--format2", dest="format2", action="store", type="string", help="format of previous file [compulsory] [format: transcript file format]")
209 parser.add_option("-5", "--upstream", dest="upstream", action="store_true", default=False, help="output upstream elements [format: boolean] [default: False]") 207 parser.add_option("-5", "--upstream", dest="upstream", action="store_true", default=False, help="output upstream elements [format: boolean] [default: False]")
210 parser.add_option("-3", "--downstream", dest="downstream", action="store_true", default=False, help="output downstream elements [format: boolean] [default: False]") 208 parser.add_option("-3", "--downstream", dest="downstream", action="store_true", default=False, help="output downstream elements [format: boolean] [default: False]")
211 parser.add_option("-c", "--colinear", dest="colinear", action="store_true", default=False, help="find first colinear element [format: boolean] [default: False]") 209 parser.add_option("-c", "--colinear", dest="colinear", action="store_true", default=False, help="find first colinear element [format: boolean] [default: False]")
212 parser.add_option("-a", "--antisense", dest="antisense", action="store_true", default=False, help="find first anti-sense element [format: boolean] [default: False]") 210 parser.add_option("-a", "--antisense", dest="antisense", action="store_true", default=False, help="find first anti-sense element [format: boolean] [default: False]")
213 parser.add_option("-e", "--noOverlap", dest="noOverlap", action="store_true", default=False, help="do not consider elements which are overlapping reference elements [format: boolean] [default: False]") 211 parser.add_option("-e", "--noOverlap", dest="noOverlap", action="store_true", default=False, help="do not consider elements which are overlapping reference elements [format: boolean] [default: False]")
214 parser.add_option("-d", "--minDistance", dest="minDistance", action="store", default=None, type="int", help="minimum distance between 2 elements [format: int]") 212 parser.add_option("-d", "--minDistance", dest="minDistance", action="store", default=None, type="int", help="minimum distance between 2 elements [format: int]")
215 parser.add_option("-D", "--maxDistance", dest="maxDistance", action="store", default=None, type="int", help="maximum distance between 2 elements [format: int]") 213 parser.add_option("-D", "--maxDistance", dest="maxDistance", action="store", default=None, type="int", help="maximum distance between 2 elements [format: int]")
216 parser.add_option("-t", "--tag", dest="tagName", action="store", default="flanking", type="string", help="name of the new tag [format: string] [default: flanking]") 214 parser.add_option("-t", "--tag", dest="tagName", action="store", default="flanking", type="string", help="name of the new tag [format: string] [default: flanking]")
217 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [format: output file in GFF3 format]") 215 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [format: output file in GFF3 format]")
218 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]") 216 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]")
219 (options, args) = parser.parse_args() 217 (options, args) = parser.parse_args()
220 218
221 gf = GetFlanking(options.verbosity) 219 gf = GetFlanking(options.verbosity)
222 gf.setInputFile(options.inputFileName1, options.format1, QUERY) 220 gf.setInputFile(options.inputFileName1, options.format1, QUERY)
223 gf.setInputFile(options.inputFileName2, options.format2, REFERENCE) 221 gf.setInputFile(options.inputFileName2, options.format2, REFERENCE)
224 gf.setOutputFile(options.outputFileName) 222 gf.setOutputFile(options.outputFileName)
225 gf.addUpstreamDirection(options.upstream) 223 gf.addUpstreamDirection(options.upstream)
226 gf.addDownstreamDirection(options.downstream) 224 gf.addDownstreamDirection(options.downstream)
227 gf.setColinear(options.colinear) 225 gf.setColinear(options.colinear)
228 gf.setAntisense(options.antisense) 226 gf.setAntisense(options.antisense)
229 gf.setNoOverlap(options.noOverlap) 227 gf.setNoOverlap(options.noOverlap)
230 gf.setMinDistance(options.minDistance) 228 gf.setMinDistance(options.minDistance)
231 gf.setMaxDistance(options.maxDistance) 229 gf.setMaxDistance(options.maxDistance)
232 gf.setNewTagName(options.tagName) 230 gf.setNewTagName(options.tagName)
233 gf.run() 231 gf.run()