|
6
|
1 #! /usr/bin/env python
|
|
|
2 #
|
|
|
3 # Copyright INRA-URGI 2009-2010
|
|
|
4 #
|
|
|
5 # This software is governed by the CeCILL license under French law and
|
|
|
6 # abiding by the rules of distribution of free software. You can use,
|
|
|
7 # modify and/ or redistribute the software under the terms of the CeCILL
|
|
|
8 # license as circulated by CEA, CNRS and INRIA at the following URL
|
|
|
9 # "http://www.cecill.info".
|
|
|
10 #
|
|
|
11 # As a counterpart to the access to the source code and rights to copy,
|
|
|
12 # modify and redistribute granted by the license, users are provided only
|
|
|
13 # with a limited warranty and the software's author, the holder of the
|
|
|
14 # economic rights, and the successive licensors have only limited
|
|
|
15 # liability.
|
|
|
16 #
|
|
|
17 # In this respect, the user's attention is drawn to the risks associated
|
|
|
18 # with loading, using, modifying and/or developing or reproducing the
|
|
|
19 # software by the user in light of its specific status of free software,
|
|
|
20 # that may mean that it is complicated to manipulate, and that also
|
|
|
21 # therefore means that it is reserved for developers and experienced
|
|
|
22 # professionals having in-depth computer knowledge. Users are therefore
|
|
|
23 # encouraged to load and test the software's suitability as regards their
|
|
|
24 # requirements in conditions enabling the security of their systems and/or
|
|
|
25 # data to be ensured and, more generally, to use and operate it in the
|
|
|
26 # same conditions as regards security.
|
|
|
27 #
|
|
|
28 # The fact that you are presently reading this means that you have had
|
|
|
29 # knowledge of the CeCILL license and that you accept its terms.
|
|
|
30 #
|
|
|
31 import re
|
|
|
32 from commons.core.writer.WriterChooser import WriterChooser
|
|
|
33 """
|
|
|
34 Cluster the data into regions (defined by size and overlap with next region) and keep only highest peaks.
|
|
|
35 """
|
|
|
36
|
|
|
37 import os, os.path
|
|
|
38 from optparse import OptionParser
|
|
|
39 from SMART.Java.Python.structure.Transcript import Transcript
|
|
|
40 from SMART.Java.Python.structure.TranscriptContainer import TranscriptContainer
|
|
|
41 from SMART.Java.Python.misc.RPlotter import RPlotter
|
|
|
42 from SMART.Java.Python.misc.Progress import Progress
|
|
|
43 from commons.core.writer.Gff3Writer import Gff3Writer
|
|
|
44
|
|
|
45 class ClusterizeBySlidingWindows(object):
|
|
|
46
|
|
|
47 def __init__(self, verbosity = 0):
|
|
|
48 self.verbosity = verbosity
|
|
|
49 self.strands = (0, )
|
|
|
50 self.normalize = False
|
|
|
51 self.plot = None
|
|
|
52 self.excel = None
|
|
|
53 self.outputFileName = ''
|
|
|
54 self.defaultValue = None
|
|
|
55
|
|
|
56 def __del__(self):
|
|
|
57 pass
|
|
|
58
|
|
|
59 def setInputFile(self, fileName, format):
|
|
|
60 self.parser = TranscriptContainer(fileName, format, self.verbosity)
|
|
|
61
|
|
|
62 def setOutputFileName(self, fileName, format="gff", title="S-MART", feature="transcript", featurePart="exon"):
|
|
|
63 writerChooser = WriterChooser(self.verbosity)
|
|
|
64 writerChooser.findFormat(format)
|
|
|
65 self.writer = writerChooser.getWriter(fileName)
|
|
|
66 self.writer.setTitle(title)
|
|
|
67 self.writer.setFeature(feature)
|
|
|
68 self.writer.setFeaturePart(featurePart)
|
|
|
69 # self.outputFileName = fileName
|
|
|
70 # self.outputFormat = format
|
|
|
71
|
|
|
72 def setWindowSize(self, size):
|
|
|
73 self.size = size
|
|
|
74
|
|
|
75 def setWindowOverlap(self, overlap):
|
|
|
76 self.overlap = overlap
|
|
|
77
|
|
|
78 def setTag(self, tag):
|
|
|
79 self.tag = tag
|
|
|
80
|
|
|
81 def setOperation(self, operation):
|
|
|
82 self.operation = operation
|
|
|
83
|
|
|
84 def setBothStrands(self, bothStrands):
|
|
|
85 if bothStrands:
|
|
|
86 self.strands = (-1, 1)
|
|
|
87
|
|
|
88 def setNormalize(self, normalize):
|
|
|
89 self.normalize = normalize
|
|
|
90
|
|
|
91 def setPlot(self, plot):
|
|
|
92 self.plot = plot
|
|
|
93
|
|
|
94 def setExcel(self, excel):
|
|
|
95 self.excel = excel
|
|
|
96
|
|
|
97 def setOutputTag(self, tag):
|
|
|
98 self.outputTagName = tag
|
|
|
99
|
|
|
100 def setDefaultValue(self, defaultValue):
|
|
|
101 self.defaultValue = defaultValue
|
|
|
102
|
|
|
103 def checkOptions(self):
|
|
|
104 # if self.operation != None:
|
|
|
105 # raise Exception("Trying to combine the values without specifying tag! Aborting...")
|
|
|
106 if self.operation != None and self.operation not in ("sum", "avg", "med", "min", "max"):
|
|
|
107 raise Exception("Do not understand tag '%s'! Aborting..." % (self.operation))
|
|
|
108
|
|
|
109 def getChromosomeSizes(self):
|
|
|
110 self.sizes = {}
|
|
|
111 progress = Progress(self.parser.getNbTranscripts(), "Getting sizes in genome", self.verbosity)
|
|
|
112 for transcript in self.parser.getIterator():
|
|
|
113 self.sizes[transcript.getChromosome()] = max(transcript.getStart(), self.sizes.get(transcript.getChromosome(), 0))
|
|
|
114 progress.inc()
|
|
|
115 progress.done()
|
|
|
116
|
|
|
117 def getBinsFromPos(self, pos):
|
|
|
118 bin = (pos - 1) / (self.size - self.overlap)
|
|
|
119 if bin >= 1 and pos <= bin * (self.size - self.overlap) + self.overlap:
|
|
|
120 return (bin - 1, bin)
|
|
|
121 return (bin, )
|
|
|
122
|
|
|
123 def getPosFromBin(self, bin):
|
|
|
124 return (bin * (self.size - self.overlap) + 1, bin * (self.size - self.overlap) + self.size)
|
|
|
125
|
|
|
126 def initializeBins(self):
|
|
|
127 self.binsPerStrand = {}
|
|
|
128 self.sumsPerStrand = {}
|
|
|
129 self.valuesPerStrand = {}
|
|
|
130 self.toBePlottedPerStrand = {}
|
|
|
131 for strand in self.strands:
|
|
|
132 self.binsPerStrand[strand] = {}
|
|
|
133 self.sumsPerStrand[strand] = {}
|
|
|
134 self.valuesPerStrand[strand] = {}
|
|
|
135 self.toBePlottedPerStrand[strand] = {}
|
|
|
136 for chromosome in self.sizes:
|
|
|
137 binRange = range(self.getBinsFromPos(self.sizes[chromosome])[-1] + 1)
|
|
|
138 self.binsPerStrand[strand][chromosome] = dict([[i, 0] for i in binRange])
|
|
|
139 self.sumsPerStrand[strand][chromosome] = dict([[i, 0.0] for i in binRange])
|
|
|
140 self.valuesPerStrand[strand][chromosome] = dict([[i, []] for i in binRange])
|
|
|
141 self.toBePlottedPerStrand[strand][chromosome] = dict([[i, 0] for i in binRange])
|
|
|
142
|
|
|
143 def getNbElements(self, transcript):
|
|
|
144 nbOccurrences = 1 if "nbOccurrences" not in transcript.getTagNames() else transcript.getTagValue("nbOccurrences")
|
|
|
145 nbElements = 1 if "nbElements" not in transcript.getTagNames() else transcript.getTagValue("nbElements")
|
|
|
146 nbOccurrences = float(nbOccurrences)
|
|
|
147 nbElements = float(nbElements)
|
|
|
148 nbElements /= float(nbOccurrences)
|
|
|
149 return nbElements
|
|
|
150
|
|
|
151 def setBins(self):
|
|
|
152 progress = Progress(self.parser.getNbTranscripts(), "Setting bins", self.verbosity)
|
|
|
153 for transcript in self.parser.getIterator():
|
|
|
154 nbElements = self.getNbElements(transcript)
|
|
|
155 strand = transcript.getDirection() if len(self.strands) == 2 else 0
|
|
|
156 for bin in self.getBinsFromPos(transcript.getStart()):
|
|
|
157 self.binsPerStrand[strand][transcript.getChromosome()][bin] += nbElements
|
|
|
158 if self.tag != None:
|
|
|
159 if self.tag not in transcript.getTagNames():
|
|
|
160 if self.defaultValue is None:
|
|
|
161 raise Exception("Tag %s undefined in transcript %s" % (self.tag, transcript))
|
|
|
162 value = self.defaultValue
|
|
|
163 else:
|
|
|
164 value = float(transcript.getTagValue(self.tag))
|
|
|
165 self.sumsPerStrand[strand][transcript.getChromosome()][bin] += value
|
|
|
166 self.valuesPerStrand[strand][transcript.getChromosome()][bin].append(value)
|
|
|
167 progress.inc()
|
|
|
168 progress.done()
|
|
|
169
|
|
|
170 def aggregateData(self):
|
|
|
171 if self.operation == "sum":
|
|
|
172 self.computeSumData()
|
|
|
173 elif self.operation == "avg":
|
|
|
174 self.computeAvgData()
|
|
|
175 elif self.operation == "med":
|
|
|
176 self.computeMedData()
|
|
|
177 elif self.operation == "min":
|
|
|
178 self.computeMinData()
|
|
|
179 elif self.operation == "max":
|
|
|
180 self.computeMaxData()
|
|
|
181 elif self.operation == "GCpercent":
|
|
|
182 self.computeGCPercent()
|
|
|
183 else:
|
|
|
184 self.toBePlottedPerStrand = self.binsPerStrand
|
|
|
185
|
|
|
186 def computeSumData(self):
|
|
|
187 self.toBePlottedPerStrand = self.sumsPerStrand
|
|
|
188
|
|
|
189 def computeAvgData(self):
|
|
|
190 for strand in self.strands:
|
|
|
191 for chromosome in self.binsPerStrand[strand]:
|
|
|
192 for bin in self.binsPerStrand[strand][chromosome]:
|
|
|
193 if self.binsPerStrand[strand][chromosome][bin] != 0:
|
|
|
194 self.toBePlottedPerStrand[strand][chromosome][bin] = float(self.sumsPerStrand[strand][chromosome][bin]) / self.binsPerStrand[strand][chromosome][bin]
|
|
|
195
|
|
|
196 def computeMedData(self):
|
|
|
197 for strand in self.strands:
|
|
|
198 for chromosome in self.binsPerStrand[strand]:
|
|
|
199 for bin in self.binsPerStrand[strand][chromosome]:
|
|
|
200 if self.valuesPerStrand[strand][chromosome][bin]:
|
|
|
201 self.valuesPerStrand[strand][chromosome][bin].sort()
|
|
|
202 size = len(self.valuesPerStrand[strand][chromosome][bin])
|
|
|
203 if size % 2 == 1:
|
|
|
204 self.toBePlottedPerStrand[strand][chromosome][bin] = self.valuesPerStrand[strand][chromosome][bin][(size - 1) / 2]
|
|
|
205 else:
|
|
|
206 self.toBePlottedPerStrand[strand][chromosome][bin] = (self.valuesPerStrand[strand][chromosome][bin][size / 2 - 1] + self.valuesPerStrand[strand][chromosome][bin][size / 2]) / 2.0
|
|
|
207
|
|
|
208 def computeMinData(self):
|
|
|
209 for strand in self.strands:
|
|
|
210 for chromosome in self.binsPerStrand[strand]:
|
|
|
211 for bin in self.binsPerStrand[strand][chromosome]:
|
|
|
212 if self.valuesPerStrand[strand][chromosome][bin]:
|
|
|
213 self.toBePlottedPerStrand[strand][chromosome][bin] = min(self.valuesPerStrand[strand][chromosome][bin])
|
|
|
214
|
|
|
215 def computeMaxData(self):
|
|
|
216 for strand in self.strands:
|
|
|
217 for chromosome in self.binsPerStrand[strand]:
|
|
|
218 for bin in self.binsPerStrand[strand][chromosome]:
|
|
|
219 if self.valuesPerStrand[strand][chromosome][bin]:
|
|
|
220 self.toBePlottedPerStrand[strand][chromosome][bin] = max(self.valuesPerStrand[strand][chromosome][bin])
|
|
|
221
|
|
|
222 def computeGCPercent(self):
|
|
|
223 for strand in self.strands:
|
|
|
224 for chromosome in self.binsPerStrand[strand]:
|
|
|
225 for bin in self.binsPerStrand[strand][chromosome]:
|
|
|
226 if self.valuesPerStrand[strand][chromosome][bin]:
|
|
|
227 subSequence = self.valuesPerStrand[strand][chromosome][bin]
|
|
|
228 NPercent = 100 * (subSequence.countNt("N") / float(subSequence.getSize()))
|
|
|
229 if NPercent >= 50:
|
|
|
230 currentGCpercent = "NA"
|
|
|
231 else:
|
|
|
232 currentGCpercent = subSequence.getGCpercentageInSequenceWithoutCountNInLength()
|
|
|
233
|
|
|
234 self.toBePlottedPerStrand[strand][chromosome][bin] = currentGCpercent
|
|
|
235 #TODO: see if a map method could be used for the various "compute" methods
|
|
|
236 #return currentGCpercent, NPercent
|
|
|
237
|
|
|
238 def plotData(self):
|
|
|
239 if self.plot != None:
|
|
|
240 for strand in self.strands:
|
|
|
241 adjunct = ""
|
|
|
242 if strand != 0:
|
|
|
243 adjunct = "Strand%d" % (strand)
|
|
|
244 for chromosome in self.toBePlottedPerStrand[strand]:
|
|
|
245 if len(self.toBePlottedPerStrand[strand][chromosome].keys()) > 0:
|
|
|
246 plotter = RPlotter(self.plot, self.verbosity)
|
|
|
247 plotter.setFill(0)
|
|
|
248 plotter.addLine(self.toBePlottedPerStrand[strand][chromosome], chromosome)
|
|
|
249 plotter.plot()
|
|
|
250
|
|
|
251 def writeExcel(self):
|
|
|
252 if self.excel != None:
|
|
|
253 excelFile = open(self.excel, "w")
|
|
|
254 for strand in self.strands:
|
|
|
255 maxBin = max([max(self.toBePlottedPerStrand[strand][chromosome].keys()) for chromosome in self.binsPerStrand[strand]])
|
|
|
256 for bin in range(0, maxBin + 1):
|
|
|
257 excelFile.write(",%d-%d" % self.getPosFromBin(bin))
|
|
|
258 excelFile.write("\n")
|
|
|
259 for chromosome in self.toBePlottedPerStrand[strand]:
|
|
|
260 excelFile.write("%s" % (chromosome))
|
|
|
261 for bin in self.toBePlottedPerStrand[strand][chromosome]:
|
|
|
262 excelFile.write(",%f" % (self.toBePlottedPerStrand[strand][chromosome][bin]))
|
|
|
263 excelFile.write("\n")
|
|
|
264 excelFile.close()
|
|
|
265
|
|
|
266 def printRegions(self):
|
|
|
267 cpt = 1
|
|
|
268 tagOp = "nb"
|
|
|
269 tagName = "Elements"
|
|
|
270 outputTagName = "nbElements"
|
|
|
271 if self.operation != None:
|
|
|
272 tagOp = self.operation.lower()
|
|
|
273 if self.tag != None:
|
|
|
274 tagName = self.tag.title()
|
|
|
275 if self.outputTagName != None:
|
|
|
276 outputTagName = self.outputTagName
|
|
|
277
|
|
|
278
|
|
|
279 #writer = Gff3Writer(self.outputFileName, self.verbosity)
|
|
|
280
|
|
|
281 for strand in self.strands:
|
|
|
282 for chromosome in self.toBePlottedPerStrand[strand]:
|
|
|
283 for bin in self.toBePlottedPerStrand[strand][chromosome]:
|
|
|
284 transcript = Transcript()
|
|
|
285 transcript.setName("region%d" % cpt)
|
|
|
286 transcript.setChromosome(chromosome)
|
|
|
287 transcript.setStart(self.getPosFromBin(bin)[0])
|
|
|
288 transcript.setEnd(self.getPosFromBin(bin)[1])
|
|
|
289 transcript.setDirection(1 if strand == 0 else strand)
|
|
|
290 transcript.setTagValue(outputTagName, self.binsPerStrand[strand][chromosome][bin])
|
|
|
291 transcript.setTagValue("%s%s" % (tagOp, tagName), str(self.toBePlottedPerStrand[strand][chromosome][bin]))
|
|
|
292 self.writer.addTranscript(transcript)
|
|
|
293 cpt += 1
|
|
|
294 self.writer.close()
|
|
|
295
|
|
|
296 def run(self):
|
|
|
297 self.checkOptions()
|
|
|
298 self.getChromosomeSizes()
|
|
|
299 self.initializeBins()
|
|
|
300 self.setBins()
|
|
|
301 self.aggregateData()
|
|
|
302 if self.excel:
|
|
|
303 self.writeExcel()
|
|
|
304 if self.plot:
|
|
|
305 self.plotData()
|
|
|
306 self.printRegions()
|
|
|
307
|
|
|
308
|
|
|
309 if __name__ == "__main__":
|
|
|
310
|
|
|
311 # parse command line
|
|
|
312 description = "Clusterize by Sliding Windows v1.0.1: Produces a GFF3 file that clusters a list of transcripts using a sliding window. [Category: Sliding Windows]"
|
|
|
313
|
|
|
314 parser = OptionParser(description = description)
|
|
|
315 parser.add_option("-i", "--input", dest="inputFileName", action="store", type="string", help="input file [compulsory] [format: file in transcript format given by -f]")
|
|
|
316 parser.add_option("-f", "--inputFormat", dest="inputFormat", action="store", type="string", help="format of the input file [compulsory] [format: transcript file format]")
|
|
|
317 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [compulsory] [format: output file in transcript format given by -u]")
|
|
|
318 parser.add_option("-u", "--outputFormat", dest="outputFormat", action="store", default="gff", type="string", help="format of the output file [format: transcript file format]")
|
|
|
319 parser.add_option("-s", "--size", dest="size", action="store", type="int", help="size of the regions [compulsory] [format: int]")
|
|
|
320 parser.add_option("-e", "--overlap", dest="overlap", action="store", type="int", help="overlap between two consecutive regions [compulsory] [format: int]")
|
|
|
321 parser.add_option("-m", "--normalize", dest="normalize", action="store_true", default=False, help="normalize the number of reads per cluster by the number of mappings per read [format: bool] [default: false]")
|
|
|
322 parser.add_option("-g", "--tag", dest="tag", action="store", default=None, type="string", help="use a given tag as input (instead of summing number of features) [format: string]")
|
|
|
323 parser.add_option("-r", "--operation", dest="operation", action="store", default=None, type="string", help="combine tag value with given operation [format: choice (sum, avg, med, min, max)]")
|
|
|
324 parser.add_option("-d", "--defaultValue",dest="defaultValue", action="store", type="float", help="default value for input tag [format: float]")
|
|
|
325 parser.add_option("-w", "--write", dest="writeTag", action="store", default=None, type="string", help="print the result in the given tag (default usually is 'nbElements') [format: string]")
|
|
|
326 parser.add_option("-2", "--strands", dest="strands", action="store_true", default=False, help="consider the two strands separately [format: bool] [default: false]")
|
|
|
327 parser.add_option("-p", "--plot", dest="plot", action="store", default=None, type="string", help="plot regions to the given file [format: output file in PNG format]")
|
|
|
328 parser.add_option("-x", "--excel", dest="excel", action="store", default=None, type="string", help="write an Excel file to the given file [format: output file in Excel format]")
|
|
|
329 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int] [default: 1]")
|
|
|
330 (options, args) = parser.parse_args()
|
|
|
331
|
|
|
332 cbsw = ClusterizeBySlidingWindows(options.verbosity)
|
|
|
333 cbsw.setInputFile(options.inputFileName, options.inputFormat)
|
|
|
334 cbsw.setOutputFileName(options.outputFileName, options.outputFormat)
|
|
|
335 cbsw.setWindowSize(options.size)
|
|
|
336 cbsw.setWindowOverlap(options.overlap)
|
|
|
337 cbsw.setTag(options.tag)
|
|
|
338 cbsw.setDefaultValue(options.defaultValue)
|
|
|
339 cbsw.setOperation(options.operation)
|
|
|
340 cbsw.setOutputTag(options.writeTag)
|
|
|
341 cbsw.setBothStrands(options.strands)
|
|
|
342 cbsw.setPlot(options.plot)
|
|
|
343 cbsw.setExcel(options.excel)
|
|
|
344 cbsw.run()
|
