Mercurial > repos > yufei-luo > s_mart
diff SMART/Java/Python/ncList/FindOverlapsWithOneInterval.py @ 36:44d5973c188c
Uploaded
author | m-zytnicki |
---|---|
date | Tue, 30 Apr 2013 15:02:29 -0400 |
parents | 769e306b7933 |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SMART/Java/Python/ncList/FindOverlapsWithOneInterval.py Tue Apr 30 15:02:29 2013 -0400 @@ -0,0 +1,197 @@ +#! /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. +# + +import struct +import math +import os +from optparse import OptionParser +from commons.core.writer.Gff3Writer import Gff3Writer +from SMART.Java.Python.ncList.NCList import NCList +from SMART.Java.Python.ncList.FileSorter import FileSorter +from commons.core.parsing.ParserChooser import ParserChooser +from SMART.Java.Python.ncList.NCListCursor import NCListCursor +from SMART.Java.Python.structure.Transcript import Transcript + +LONGSIZE = struct.calcsize('l') + +class FindOverlapsWithOneInterval(object): + + def __init__(self, verbosity): + self._sortedFileName = None + self._verbosity = verbosity + self._overlappingNames = [] + self._nbOverlaps = 0 + self._nbWritten = 0 + + def __del__(self): + if self._sortedFileName and os.path.exists(self._sortedFileName): + os.remove(self._sortedFileName) + + def close(self): + self._iWriter.close() + + def setOutputFileName(self, fileName): + self._iWriter = Gff3Writer(fileName) + + def setFileName(self, fileName, format): + chooser = ParserChooser(self._verbosity) + chooser.findFormat(format) + self._parser = chooser.getParser(fileName) + self._sortedFileName = "%s_sorted.pkl" % (os.path.splitext(fileName)[0]) + + def setInterval(self, chromosome, start, end): + self._chromosome = chromosome + self._start = start + self._end = end + self._transcript = Transcript() + self._transcript.setChromosome(chromosome) + self._transcript.setStart(start) + self._transcript.setEnd(end) + self._transcript.setDirection("+") + + def setTranscript(self, transcript): + if transcript.__class__.__name__ == "Mapping": + transcript = transcript.getTranscript() + self._chromosome = transcript.getChromosome() + self._start = transcript.getStart() + self._end = transcript.getEnd() + self._transcript = transcript + + def prepareIntermediateFiles(self): + fs = FileSorter(self._parser, self._verbosity-4) + fs.selectChromosome(self._chromosome) + fs.perChromosome(False) + fs.setOutputFileName(self._sortedFileName) + fs.sort() + self._nbTotalLines = fs.getNbElements() + self._nbLines = fs.getNbElementsPerChromosome()[self._chromosome] + + def createNCList(self): + if self._verbosity > 2: + print "Creating NC-list..." + ncList = NCList(self._verbosity) + ncList.createIndex(True) + ncList.setChromosome(self._chromosome) + ncList.setFileName(self._sortedFileName) + ncList.setNbElements(self._nbTotalLines) + ncList.buildLists() + self.setNCList(ncList, ncList.getIndex()) + if self._verbosity > 2: + print " ...done (%ds)" % (endTime - startTime) + + def setNCList(self, ncList, index): + self._ncList = ncList + self._indix = index + + def binarySearch(self, cursor, startL, endL): + if startL > endL: + return None + middleL = (startL + endL) / 2 + cursor.moveSibling(middleL) + overlap = self.isOverlapping(cursor) + if overlap == 0: + if middleL == startL: + return cursor + else: + return self.binarySearch(cursor, startL, middleL) + if overlap == -1: + return self.binarySearch(cursor, middleL + 1, endL) + return self.binarySearch(cursor, startL, middleL - 1) + + def compare(self, cursor = None): + self._ncList.openFiles() + if cursor == None: + dump = True + cursor = NCListCursor(None, self._ncList, 0, self._verbosity) + cursor._getSiblingData() + cursor = self.binarySearch(cursor, cursor._firstSiblingLIndex, cursor._lastSiblingLIndex) + if cursor == None: + return + while not cursor.isOut() and self.isOverlapping(cursor) == 0: + self.write(cursor) + newCursor = NCListCursor(cursor) + if newCursor.hasChildren(): + newCursor.moveDown() + self.compare(newCursor) + if cursor.isLast(): + return + cursor.moveRight() + + def isOverlapping(self, cursor): + if self._end < cursor.getStart(): + return 1 + if self._start > cursor.getEnd(): + return -1 + return 0 + + def write(self, cursor): + self._nbOverlaps += 1 + refTranscript = cursor.getTranscript() + self._overlappingNames.append(refTranscript.getName()) + + def dumpWriter(self): + if (not self._overlappingNames) or self._transcript == None: + return + self._transcript.setTagValue("nbOverlaps", len(self._overlappingNames)) + self._transcript.setTagValue("overlapsWith", "--".join(self._overlappingNames)) + self._iWriter.addTranscript(self._transcript) + self._nbWritten += 1 + self._overlappingNames = [] + + def run(self): + self.prepareIntermediateFiles() + self.createNCList() + self.compare() + self.dumpWriter() + self.close() + if self._verbosity > 0: + print "# refs: %d" % (self._nbLines) + print "# written: %d (%d overlaps)" % (self._nbOverlappingQueries, self._nbOverlaps) + + +if __name__ == "__main__": + description = "FindOverlapsWithOneInterval: Finds overlaps with one query interval." + + 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 previous file [compulsory] [format: transcript file format]") + parser.add_option("-s", "--start", dest="start", action="store", type="int", help="The start of the query interval [compulsory] [format: int]") + parser.add_option("-e", "--end", dest="end", action="store", type="int", help="The end of the query interval [compulsory] [format: int]") + parser.add_option("-c", "--chromosome", dest="chromosome", action="store", type="string", help="Chromosome of the query interval [compulsory] [format: string]") + parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="Output file [compulsory] [format: output file in GFF3 format]") + 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() + + iFOWOI = FindOverlapsWithOneInterval(options.verbosity) + iFOWOI.setFileName(options.inputFileName, options.format) + iFOWOI.setInterval(options.chromosome, options.start, options.end) + iFOWOI.setOutputFileName(options.outputFileName) + iFOWOI.run()