Mercurial > repos > yufei-luo > s_mart
diff smart_toolShed/SMART/Java/Python/clusterize.py @ 0:e0f8dcca02ed
Uploaded S-MART tool. A toolbox manages RNA-Seq and ChIP-Seq data.
author | yufei-luo |
---|---|
date | Thu, 17 Jan 2013 10:52:14 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/smart_toolShed/SMART/Java/Python/clusterize.py Thu Jan 17 10:52:14 2013 -0500 @@ -0,0 +1,165 @@ +#! /usr/bin/env python +# +# Copyright INRA-URGI 2009-2010 +# +# This software is governed by the CeCILL license under French law and +# abiding by the rules of distribution of free software. You can use, +# modify and/ or redistribute the software under the terms of the CeCILL +# license as circulated by CEA, CNRS and INRIA at the following URL +# "http://www.cecill.info". +# +# As a counterpart to the access to the source code and rights to copy, +# modify and redistribute granted by the license, users are provided only +# with a limited warranty and the software's author, the holder of the +# economic rights, and the successive licensors have only limited +# liability. +# +# In this respect, the user's attention is drawn to the risks associated +# with loading, using, modifying and/or developing or reproducing the +# software by the user in light of its specific status of free software, +# that may mean that it is complicated to manipulate, and that also +# therefore means that it is reserved for developers and experienced +# professionals having in-depth computer knowledge. Users are therefore +# encouraged to load and test the software's suitability as regards their +# requirements in conditions enabling the security of their systems and/or +# data to be ensured and, more generally, to use and operate it in the +# same conditions as regards security. +# +# The fact that you are presently reading this means that you have had +# knowledge of the CeCILL license and that you accept its terms. +# +from commons.core.writer.WriterChooser import WriterChooser +"""Clusterize a set of transcripts""" + +import os +from optparse import OptionParser +from commons.core.parsing.ParserChooser import ParserChooser +from commons.core.writer.Gff3Writer import Gff3Writer +from SMART.Java.Python.structure.Transcript import Transcript +from SMART.Java.Python.ncList.NCListFilePickle import NCListFileUnpickle +from SMART.Java.Python.ncList.FileSorter import FileSorter +from SMART.Java.Python.misc.Progress import Progress + +class Clusterize(object): + + def __init__(self, verbosity): + self.normalize = False + self.presorted = False + self.distance = 1 + self.colinear = False + self.nbWritten = 0 + self.nbMerges = 0 + self.verbosity = verbosity + self.splittedFileNames = {} + + def __del__(self): + for fileName in self.splittedFileNames.values(): + os.remove(fileName) + + def setInputFile(self, fileName, format): + parserChooser = ParserChooser(self.verbosity) + parserChooser.findFormat(format) + self.parser = parserChooser.getParser(fileName) + self.sortedFileName = "%s_sorted.pkl" % (os.path.splitext(fileName)[0]) + + def setOutputFileName(self, fileName, format="gff3", title="S-MART", feature="transcript", featurePart="exon"): + writerChooser = WriterChooser() + writerChooser.findFormat(format) + self.writer = writerChooser.getWriter(fileName) + self.writer.setTitle(title) + self.writer.setFeature(feature) + self.writer.setFeaturePart(featurePart) + + def setDistance(self, distance): + self.distance = distance + + def setColinear(self, colinear): + self.colinear = colinear + + def setNormalize(self, normalize): + self.normalize = normalize + + def setPresorted(self, presorted): + self.presorted = presorted + + def _sortFile(self): + fs = FileSorter(self.parser, self.verbosity-4) + fs.perChromosome(True) + fs.setPresorted(self.presorted) + fs.setOutputFileName(self.sortedFileName) + fs.sort() + self.splittedFileNames = fs.getOutputFileNames() + self.nbElementsPerChromosome = fs.getNbElementsPerChromosome() + self.nbElements = fs.getNbElements() + + def _iterate(self, chromosome): + progress = Progress(self.nbElementsPerChromosome[chromosome], "Checking chromosome %s" % (chromosome), self.verbosity) + transcripts = [] + parser = NCListFileUnpickle(self.splittedFileNames[chromosome], self.verbosity) + for newTranscript in parser.getIterator(): + newTranscripts = [] + for oldTranscript in transcripts: + if self._checkOverlap(newTranscript, oldTranscript): + self._merge(newTranscript, oldTranscript) + elif self._checkPassed(newTranscript, oldTranscript): + self._write(oldTranscript) + else: + newTranscripts.append(oldTranscript) + newTranscripts.append(newTranscript) + transcripts = newTranscripts + progress.inc() + for transcript in transcripts: + self._write(transcript) + progress.done() + + def _merge(self, transcript1, transcript2): + self.nbMerges += 1 + transcript2.setDirection(transcript1.getDirection()) + transcript1.merge(transcript2) + + def _write(self, transcript): + self.nbWritten += 1 + self.writer.addTranscript(transcript) + + def _checkOverlap(self, transcript1, transcript2): + if self.colinear and transcript1.getDirection() != transcript2.getDirection(): + return False + if transcript1.getDistance(transcript2) > self.distance: + return False + return True + + def _checkPassed(self, transcript1, transcript2): + return (transcript1.getDistance(transcript2) > self.distance) + + def run(self): + self._sortFile() + for chromosome in sorted(self.splittedFileNames.keys()): + self._iterate(chromosome) + self.writer.close() + if self.verbosity > 0: + print "# input: %d" % (self.nbElements) + print "# written: %d (%d%% overlaps)" % (self.nbWritten, 0 if (self.nbElements == 0) else ((float(self.nbWritten) / self.nbElements) * 100)) + print "# merges: %d" % (self.nbMerges) + + +if __name__ == "__main__": + description = "Clusterize v1.0.3: clusterize the data which overlap. [Category: Merge]" + + parser = OptionParser(description = description) + parser.add_option("-i", "--input", dest="inputFileName", action="store", type="string", help="input file [compulsory] [format: file in transcript format given by -f]") + parser.add_option("-f", "--format", dest="format", action="store", type="string", help="format of file [format: transcript file format]") + parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [compulsory] [format: output file in transcript format given by -u]") + parser.add_option("-u", "--outputFormat", dest="outputFormat", action="store", default="gff", type="string", help="output file format [format: transcript file format]") + parser.add_option("-c", "--colinear", dest="colinear", action="store_true", default=False, help="merge colinear transcripts only [format: bool] [default: false]") + parser.add_option("-d", "--distance", dest="distance", action="store", default=0, type="int", help="max. distance between two transcripts to be merged [format: int] [default: 0]") + parser.add_option("-n", "--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]") + parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int] [default: 1]") + (options, args) = parser.parse_args() + + c = Clusterize(options.verbosity) + c.setInputFile(options.inputFileName, options.format) + c.setOutputFileName(options.outputFileName, options.outputFormat) + c.setColinear(options.colinear) + c.setDistance(options.distance) + c.setNormalize(options.normalize) + c.run()