Mercurial > repos > yufei-luo > s_mart
diff SMART/Java/Python/structure/Interval.py @ 36:44d5973c188c
Uploaded
author | m-zytnicki |
---|---|
date | Tue, 30 Apr 2013 15:02:29 -0400 |
parents | |
children | 169d364ddd91 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SMART/Java/Python/structure/Interval.py Tue Apr 30 15:02:29 2013 -0400 @@ -0,0 +1,706 @@ +# +# 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 SMART.Java.Python.structure.Bins import * +from commons.core.coord.Range import Range + +class Interval(Range): + """ + Store a genomic interval + @ivar name: name of the interval [optional] + @type name: string + @ivar id: id of the interval [optional] + @type id: int + @ivar bin: bin in which the interval should be if stored in a database [computed] + @type bin: int + @ival tags: information about the transcript [optional] + @type tags: dict + @ivar verbosity: verbosity + @type verbosity: int [default: 0] + """ + + def __init__(self, interval = None, verbosity = 0): + """ + Constructor + @param interval: interval to be copied + @type interval: class L{Interval<Interval>} + @param verbosity: verbosity + @type verbosity: int + """ + Range.__init__(self) + self.name = None + self.id = None + self.bin = None + self.verbosity = verbosity + self.tags = {} + if interval != None: + self.copy(interval) + + #!!!! Warning: two methods getStart() and getEnd() give the information maximum and minimum in interval.!!!!# + #In case strand = "+", start < end; strand = "-", start > end + def getStart(self): + if self.start == -1: + return -1 + if self.end == -1: + return self.start + return self.getMin() + + + def getEnd(self): + if self.end == -1: + return -1 + if self.start == -1: + return self.end + return self.getMax() + + + def getChromosome(self): + return self.getSeqname() + + + def getDirection(self): + return 1 if self.getStrand() == "+" else -1 + + + def getName(self): + return self.name + + + def isSet(self): + """ + Check if the interval is set + """ + return self.getStart() == None and self.getEnd() == None + + + def copy(self, interval): + """ + Copy method + @param interval: interval to be copied + @type interval: class L{Interval<Interval>} + """ + self.setStart(interval.getStart()) + self.setEnd(interval.getEnd()) + self.setChromosome(interval.getChromosome()) + self.setDirection(interval.getDirection()) + self.name = interval.name + self.id = interval.id + self.bin = interval.bin + self.tags = {} + for tag in interval.tags: + self.tags[tag] = interval.tags[tag] + self.verbosity = interval.verbosity + + + def setName(self, name): + """ + Set the name + @param name: name of the interval + @type name: string + """ + if len(name) > 100: + name = name[:100] + self.name = name + + + def setChromosome(self, chromosome=""): + """ + Set the chromosome + @param chromosome: chromosome on which the interval is + @type chromosome: string + """ + if not chromosome: + self.seqname = None + else: + self.seqname = chromosome.replace(".", "_").replace("|", "_") + + + def setStart(self, start): + """ + Set the start point + Possibly reset bin + @param start: start point of the interval + @type start: int + """ + self.bin = None + direction = self.getDirection() + if self.start == -1: + self.start = start + elif self.end == -1: + self.end = start + else: + if direction == 1: + self.start = start + else: + self.end = start + if direction == 1: + self.start, self.end = min(self.start, self.end), max(self.start, self.end) + else: + self.start, self.end = max(self.start, self.end), min(self.start, self.end) + + + def setEnd(self, end): + """ + Set the end point + Possibly reset bin + @param end: end point of the interval of the interval + @type end: int + """ + self.bin = None + direction = self.getDirection() + if self.end == -1: + self.end = end + elif self.start == -1: + self.start = end + else: + if direction == 1: + self.end = end + else: + self.start = end + if direction == 1: + self.start, self.end = min(self.start, self.end), max(self.start, self.end) + else: + self.start, self.end = max(self.start, self.end), min(self.start, self.end) + + + def setSize(self, size): + """ + Possibly modify the end point + @param size: size of the transcript + @type size: int + """ + if self.end == None and self.start != None: + self.setEnd(self.start + self.getSize() - 1) + elif self.start == None and self.end != None: + self.setStart(self.end - self.getSize() + 1) + + + def getSize(self): + """ + Get the size + """ + return self.getEnd() - self.getStart() + 1 + + + def _setDirection(self, direction): + """ + Set the direction of the interval (connection to Range) + @param direction: direction of the transcript (+ / -) + @type direction: int (1 or -1) + """ + if direction * self.getDirection() < 0: + self.reverse() + + + def setDirection(self, direction): + """ + Set the direction of the interval + Possibly parse different formats + @param direction: direction of the transcript (+ / -) + @type direction: int or string + """ + if type(direction).__name__ == 'int': + self._setDirection(direction / abs(direction)) + elif type(direction).__name__ == 'str': + if direction == "+": + self._setDirection(1) + elif direction == "-": + self._setDirection(-1) + elif direction == "1" or direction == "-1": + self._setDirection(int(direction)) + elif direction.lower() == "plus": + self._setDirection(1) + elif direction.lower() == "minus": + self._setDirection(-1) + else: + raise Exception("Cannot understand direction %s" % (direction)) + else: + raise Exception("Cannot understand direction %s" % (direction)) + + + def extendStart(self, size): + """ + Extend the interval by the 5' end + @param size: the size to be exended + @type size: int + """ + if self.getDirection() == 1: + self.setStart(max(0, self.getStart() - size)) + else: + self.setEnd(self.getEnd() + size) + self.bin = None + + + def extendEnd(self, size): + """ + Extend the interval by the 3' end + @param size: the size to be exended + @type size: int + """ + if self.getDirection() == 1: + self.setEnd(self.getEnd() + size) + else: + self.setStart(max(0, self.getStart() - size)) + self.bin = None + + + def restrictStart(self, size = 1): + """ + Restrict the interval by some nucleotides, start from its start position + Remove the exons + @param size: the size to be restricted to + @type size: int + """ + if self.getDirection() == 1: + self.setEnd(min(self.getEnd(), self.getStart() + size - 1)) + else: + self.setStart(max(self.getStart(), self.getEnd() - size + 1)) + self.bin = None + + + def restrictEnd(self, size = 1): + """ + Restrict the interval by some nucleotides, end from its end position + Remove the exons + @param size: the size to be restricted to + @type size: int + """ + if self.getDirection() == 1: + self.setStart(max(self.getStart(), self.getEnd() - size + 1)) + else: + self.setEnd(min(self.getEnd(), self.getStart() + size - 1)) + self.bin = None + + + + def setTagValue(self, name, value): + """ + Set a tag + @param name: name of the tag + @type name: string + @param value: value of the tag + @type value: int or string + """ + self.tags[name] = value + + + def getTagNames(self): + """ + Get all the names of the tags + """ + return self.tags.keys() + + + def getTagValue(self, tag): + """ + Get the value of a tag + @param tag: name of a tag + @type tag: string + """ + if tag not in self.tags: + return None + return self.tags[tag] + + + def getTagValues(self, tagSep = "; ", fieldSep = " ", surrounder = ""): + """ + Get the formatted tag values + @param tagSep: separator between tags + @type tagSep: string + @param fieldSep: separator between tag name and tag value + @type fieldSep: string + @param surrounder: string which optionally surround values + @type surrounder: string + """ + tags = [] + for name, value in self.tags.iteritems(): + if value == None: + continue + if isinstance(value, basestring): + tags.append("%s%s%s%s%s" % (name, fieldSep, surrounder, value.replace("'", "\\'"), surrounder)) + elif type(value) is int: + tags.append("%s%s%s%i%s" % (name, fieldSep, surrounder, value, surrounder)) + elif type(value) is float: + tags.append("%s%s%s%f%s" % (name, fieldSep, surrounder, value, surrounder)) + else: + raise Exception("Do not know how to print '" + value + "'.") + if self.getName() != None: + tags.append("%s%s%s%s%s" % ("Name", fieldSep, surrounder, self.getName(), surrounder)) + return tagSep.join(tags) + + + def setTagValues(self, tags, tagSep = "; ", fieldSep = " "): + """ + Set the tag values using given string + @param tags: the tags, concatenated + @type tags: string + @param tagSep: separator between tags + @type tagSep: string + @param fieldSep: separator between tag name and tag value + @type fieldSep: string + """ + if tags == "": + self.tags = {} + return + for splittedTag in tags.split(tagSep): + if fieldSep not in splittedTag: + raise Exception("Weird field '%s' in tags '%s'" % (splittedTag, tags)) + tag, value = splittedTag.split(fieldSep, 1) + if tag == "Name": + self.setName(value) + continue + try: + intValue = int(value) + self.tags[tag] = intValue + except ValueError: + try: + floatValue = float(value) + self.tags[tag] = floatValue + except ValueError: + self.tags[tag] = value + + + def deleteTag(self, tag): + """ + Remove a tag + @param tag: the tag to be removed + @type tag: string + """ + if tag in self.tags: + del self.tags[tag] + + + def setNbOccurrences(self, nbOccurrences): + """ + Set the number of occurrences of the interval + @param nbOccurrences: number of occurrences of the interval + @type nbOccurrences: int + """ + self.setTagValue("nbOccurrences", nbOccurrences) + + + def setOccurrence(self, occurrence): + """ + Set the occurrence of this interval + @param occurrence: an occurrence for this transcript + @type occurrence: int + """ + self.setTagValue("occurrence", occurrence) + + def __eq__(self, interval): + """ + Whether two intervals are equal (start and end at same position) + @param interval: object to be compared to + @type interval: class L{Interval<Interval>} + """ + if not interval: + return False + return self.getChromosome() == interval.getChromosome() and self.getStart() == interval.getStart() and self.getEnd() == interval.getEnd() and self.getDirection() == interval.getDirection() + + + def overlapWith(self, interval, nbNucleotides = 1): + """ + Whether two intervals overlap + @param interval: object to be compared to + @type interval: class L{Interval<Interval>} + @param nbNucleotides: minimum number of nucleotides to declare and overlap + @type nbNucleotides: int + """ + if self.getChromosome() != interval.getChromosome(): + return False + return (min(self.getEnd(), interval.getEnd()) - max(self.getStart(), interval.getStart()) + 1 >= nbNucleotides) + + def isIncludeIn(self, interval): + return interval.include(self) + + + def include(self, interval): + """ + Whether this interval includes the other one + @param interval: object to be compared to + @type interval: class L{Interval<Interval>} + """ + if self.getChromosome() != interval.getChromosome(): + return False + return ((self.getStart() <= interval.getStart()) and (self.getEnd() >= interval.getEnd())) + + + def getDifference(self, interval, sameStrand = False): + """ + Get the difference between this cluster and another one + @param interval: object to be compared to + @type interval: class L{Interval<Interval>} + @param sameStrand: do the comparison iff the intervals are on the same strand + @type sameStrand: boolean + @return: a (possibly empty) list of intervals + """ + newInterval = Interval() + newInterval.copy(self) + if self.getChromosome() != interval.getChromosome(): + return [newInterval] + if not self.overlapWith(interval): + return [newInterval] + if sameStrand and self.getDirection() != interval.getDirection(): + return [newInterval] + intervals = [] + if self.getStart() < interval.getStart(): + newInterval = Interval() + newInterval.copy(self) + newInterval.setEnd(min(self.getEnd(), interval.getStart() - 1)) + intervals.append(newInterval) + if self.getEnd() > interval.getEnd(): + newInterval = Interval() + newInterval.copy(self) + newInterval.setStart(max(self.getStart(), interval.getEnd() + 1)) + intervals.append(newInterval) + return intervals + + + def getIntersection(self, interval): + """ + Get the intersection between this interval and another one + @param interval: object to be compared to + @type interval: class L{Interval<Interval>} + @return: an other interval + """ + if not self.overlapWith(interval): + return None + newInterval = Interval() + newInterval.setChromosome(self.getChromosome()) + newInterval.setDirection(self.getDirection()) + newInterval.setName("%s_intersect_%s" % (self.getName(), interval.getName())) + newInterval.setStart(max(self.getStart(), interval.getStart())) + newInterval.setEnd(min(self.getEnd(), interval.getEnd())) + return newInterval + + + def getDistance(self, interval): + """ + Get the distance between two intervals (a non-negative value) + @param interval: another interval + @type interval: class L{Interval<Interval>} + """ + if self.overlapWith(interval): + return 0 + if self.getChromosome() != interval.getChromosome(): + raise Exception("Cannot get the distance between %s and %s" % (str(self), str(interval))) + return min(abs(self.getStart() - interval.getEnd()), abs(self.getEnd() - interval.getStart())) + + + def getRelativeDistance(self, interval): + """ + Get the distance between two intervals (negative if first interval is before) + @param interval: another interval + @type interval: class L{Interval<Interval>} + """ + if self.overlapWith(interval): + return 0 + if self.getChromosome() != interval.getChromosome(): + raise Exception("Cannot get the distance between %s and %s" % (str(self), str(interval))) + if self.getEnd() < interval.getStart(): + distance = interval.getStart() - self.getEnd() + else: + distance = interval.getEnd() - self.getStart() + distance *= self.getDirection() + return distance + + + def merge(self, interval, normalization = False): + """ + Merge two intervals + @param interval: another interval + @type interval: class L{Interval<Interval>} + @param normalization: whether the sum of the merge should be normalized wrt the number of mappings of each elements + @type normalization: boolean + """ + if self.getChromosome() != interval.getChromosome(): + raise Exception("Cannot merge '%s' and '%s' for they are on different chromosomes." % (str(self), str(interval))) + direction = None + if self.getStart() == self.getEnd(): + direction = interval.getDirection() + elif interval.getStart() == interval.getEnd(): + direction = self.getDirection() + elif self.getDirection() != interval.getDirection(): + raise Exception("Cannot merge '%s' and '%s' for they are on different strands." % (str(self), str(interval))) + self.setStart(min(self.getStart(), interval.getStart())) + self.setEnd(max(self.getEnd(), interval.getEnd())) + if direction != None: + self.setDirection(direction) + nbElements = 0.0 + for element in (self, interval): + for tagName in ("nbElements", "nbOccurrences"): + if tagName not in element.getTagNames(): + element.setTagValue(tagName, 1) + nbElements += float(element.getTagValue("nbElements")) / float(element.getTagValue("nbOccurrences")) if normalization else float(element.getTagValue("nbElements")) + self.setTagValue("nbElements", nbElements) + self.bin = None + for tagName in ("identity", "nbOccurrences", "occurrence", "nbMismatches", "nbGaps", "rank", "evalue", "bestRegion"): + if tagName in self.getTagNames(): + del self.tags[tagName] + + + def getBin(self): + """ + Get the bin of the interval + Computed on the fly + """ + if self.bin == None: + self.bin = getBin(self.getStart(), self.getEnd()) + return self.bin + + + def getBins(self): + """ + Get all the bin this interval could fall into + """ + return getOverlappingBins(self.getStart(), self.getEnd()) + + + def getSqlVariables(cls): + """ + Get the properties of the object that should be saved in a database + """ + variables = ["name", "chromosome", "start", "end", "direction", "tags", "bin"] + return variables + getSqlVariables = classmethod(getSqlVariables) + + + def setSqlValues(self, array): + """ + Set the values of the properties of this object as given by a results line of a SQL query + """ + self.id = array[0] + self.name = array[1].strip("'") + self.setChromosome(array[2].strip("'")) + self.setStart(array[3]) + self.setEnd(array[4]) + self.setDirection(array[5]) + self.setTagValues(array[6].strip("'"), ";", "=") + self.bin = array[7] + + + def getSqlValues(self): + """ + Get the values of the properties that should be saved in a database + """ + values = dict() + values["name"] = self.name + values["chromosome"] = self.getChromosome() + values["start"] = self.getStart() + values["end"] = self.getEnd() + values["direction"] = self.getDirection() + values["tags"] = self.getTagValues(";", "=") + values["bin"] = self.getBin() + return values + + + def getSqlTypes(cls): + """ + Get the values of the properties that should be saved in a database + """ + types = dict() + types["name"] = "varchar" + types["chromosome"] = "varchar" + types["start"] = "int" + types["end"] = "int" + types["direction"] = "tinyint" + types["tags"] = "varchar" + types["bin"] = "int" + return types + getSqlTypes = classmethod(getSqlTypes) + + + def getSqlSizes(cls): + """ + Get the sizes of the properties that should be saved in a database + """ + sizes = dict() + sizes["name"] = 255 + sizes["chromosome"] = 255 + sizes["start"] = 11 + sizes["end"] = 11 + sizes["direction"] = 4 + sizes["tags"] = 1023 + sizes["bin"] = 11 + return sizes + getSqlSizes = classmethod(getSqlSizes) + + + def printCoordinates(self): + """ + Print start and end positions (depending on the direction of the interval) + """ + if self.getDirection() == 1: + return "%d-%d" % (self.getStart(), self.getEnd()) + else: + return "%d-%d" % (self.getEnd(), self.getStart()) + + + def extractSequence(self, parser): + """ + Get the sequence corresponding to this interval + @param parser: a parser to a FASTA file + @type parser: class L{SequenceListParser<SequenceListParser>} + @return : a instance of L{Sequence<Sequence>} + """ + return parser.getSubSequence(self.getChromosome(), self.getStart(), self.getEnd(), self.getDirection(), self.name) + + + def extractWigData(self, parser): + """ + Get the data retrieved from a wig file + @param parser: a parser class to a WIG file + @type parser: class L{WigParser<WigParser>} + """ + data = parser.getRange(self.getChromosome(), self.getStart(), self.getEnd()) + if self.getDirection() == -1: + if parser.strands: + newData = {} + for strand in data: + data[strand].reverse() + newData[-strand] = data[strand] + data = newData + else: + data.reverse() + return data + + + def __str__(self): + """ + Output a simple representation of this interval + """ + direction = "+" + if self.getDirection() == -1: + direction = "-" + string = "%s:%d-%d (%s)" % (self.getChromosome(), self.getStart(), self.getEnd(), direction) + if self.name != "": + string = "(%s) %s" % (self.name, string) + return string +