annotate SMART/Java/Python/structure/TranscriptListsComparator.py @ 53:47310c4fb725

Uploaded

1 #
2 # Copyright INRA-URGI 2009-2010
3 # 
4 # This software is governed by the CeCILL license under French law and
5 # abiding by the rules of distribution of free software. You can use,
6 # modify and/ or redistribute the software under the terms of the CeCILL
7 # license as circulated by CEA, CNRS and INRIA at the following URL
8 # "http://www.cecill.info".
9 # 
10 # As a counterpart to the access to the source code and rights to copy,
11 # modify and redistribute granted by the license, users are provided only
12 # with a limited warranty and the software's author, the holder of the
13 # economic rights, and the successive licensors have only limited
14 # liability.
15 # 
16 # In this respect, the user's attention is drawn to the risks associated
17 # with loading, using, modifying and/or developing or reproducing the
18 # software by the user in light of its specific status of free software,
19 # that may mean that it is complicated to manipulate, and that also
20 # therefore means that it is reserved for developers and experienced
21 # professionals having in-depth computer knowledge. Users are therefore
22 # encouraged to load and test the software's suitability as regards their
23 # requirements in conditions enabling the security of their systems and/or
24 # data to be ensured and, more generally, to use and operate it in the
25 # same conditions as regards security.
26 # 
27 # The fact that you are presently reading this means that you have had
28 # knowledge of the CeCILL license and that you accept its terms.
29 #
30 import sys
31 import random
32 from SMART.Java.Python.misc import Utils
33 from SMART.Java.Python.structure.Transcript import Transcript
34 from SMART.Java.Python.structure.TranscriptList import TranscriptList
35 from SMART.Java.Python.structure.TranscriptContainer import TranscriptContainer
36 from SMART.Java.Python.mySql.MySqlConnection import MySqlConnection
37 from SMART.Java.Python.mySql.MySqlTranscriptTable import MySqlTranscriptTable
38 from SMART.Java.Python.misc.Progress import Progress
39 from commons.core.writer.MySqlTranscriptWriter import MySqlTranscriptWriter
40 
41 
42 
43 class TranscriptListsComparator(object):
44     """
45     Compare two transcript lists, using a database for one of the list
46     Uses one TranscriptContainer for query data, 
47              one TranscriptContainer exported to MySqlTranscriptTable for reference data, 
48              one MySqlTranscriptTable for transformed reference data
49     @ivar inputTranscriptContainers: parsers to the list of query transcripts
50     @type inputTranscriptContainers: list of 2 L{TranscriptContainer<TranscriptContainer>}
51     @ivar writer:                    transcript list writer
52     @type writer:                    class L{TranscriptListWriter<TranscriptListWriter>}
53     @ivar mySqlConnection:           connection to a MySQL database (to compute the ovelapping efficiently)
54     @type mySqlConnection:           class L{MySqlConnection<MySqlConnection>}
55     @ivar introns:                   compare transcripts or exons only
56     @type introns:                   list of 2 boolean
57     @ivar starts:                    restrict the query transcripts to first nucleotides
58     @type starts:                    list of 2 int or None
59     @ivar fivePrimes:                extend a list of transcripts by their 5' end
60     @type fivePrimes:                list of 2 int or None
61     @ivar threePrimes:               extend a list of transcripts by their 3' end
62     @type threePrimes:               list of 2 int or None
63     @ivar minDistance:               min distance between two transcripts [default: 0]
64     @type minDistance:               int
65     @ivar maxDistance:               max distance between two transcripts [default: 0]
66     @type maxDistance:               int
67     @ivar minOverlap:                minimum number of overlapping nucleotides to declare an overlap
68     @type minOverlap:                int
69     @ivar pcOverlap:                 percentage of overlapping nucleotides to declare an overlap
70     @type pcOverlap:                 int
71     @ivar upstreams:                 consider distances with elements which are upstream of the transcripts
72     @type upstreams:                 boolean
73     @ivar downstreams:               consider distances with elements which are downstream of the transcripts
74     @type downstreams:               boolean
75     @ivar colinear:                  whether transcripts should overlap in the same direction
76     @type colinear:                  boolean
77     @ivar antisense:                 whether transcripts should overlap in the opposite direction
78     @type antisense:                 boolean
79     @ivar outputDistance:            output distance between query and reference instead of query transcript
80     @type outputDistance:            boolean
81     @ivar absolute:                  do not consider the strand while computing distance
82     @type absolute:                  boolean
83     @ivar strandedDistance:          return a line per strand while computing distances
84     @type strandedDistance:          boolean
85     @ivar QUERY:                     constant specifying the query objects
86     @type QUERY:                     int
87     @ivar REFERENCE:                 constant specifying the reference objects
88     @type REFERENCE:                 int
89     @ivar INPUTTYPES:                set of input types of data (query or reference) objects
90     @type INPUTTYPES:                list of 2 int
91     @ivar typeToString:              string representation of the previous types
92     @type typeToString:              dict
93     @ivar tableNames:                name of the transcript tables
94     @type tableNames:                dict of strings
95     @ivar nbTranscripts:             number of transcript in the query/reference set
96     @type nbTranscripts:             list of 2 int or None
97     @ivar nbNucleotides:             number of nucleotides in the query/reference set
98     @type nbNucleotides:             list of 2 int or None
99     @ivar transcriptsToBeStored:     transcripts that will be stored into database
00     @type transcriptsToBeStored:     dict of class L{TranscriptList<TranscriptList>}
01     @ivar multiple:                  in merge mode, aggregate multiple transcripts
02     @type multiple:                  boolean
03     @ivar normalization:             normalize each element by the number of mappings of this element
04     @type normalization:             boolean
05     @ivar invert:                    invert the current comparison
06     @type invert:                    boolean
07     @ivar splitDifference:           split into intervals when computing difference
08     @type splitDifference:           boolean
09     @ivar odds:                      whether odds about the comparison should be computed
10     @type odds:                      boolean
11     @ivar overlapResults:            count the number of overlaps
12     @type overlapResults:            dictionary
13     @ivar oddResults:                compute the number of times each transcript overlaps (or is merged with) another one
14     @type oddResults:                dictionary
15     @ivar outputContainer:           container of the output transcripts
16     @type outputContainer:           class L{TranscriptContainer<TranscriptContainer>}
17     @ivar logHandle:                 log handle
18     @type logHandle:                 file
19     @ivar verbosity:                 verbosity
20     @type verbosity:                 int    
21     """
22     
23     def __init__(self, logHandle = None, verbosity = 0):
24         """
25         Constructor
26         @param transcriptListParser2: parser to the list of reference transcripts
27         @type  transcriptListParser2: class L{TranscriptListParser<TranscriptListParser>}
28         @param logHandle:             log handle
29         @type  logHandle:             file
30         @param verbosity:             verbosity
31         @type  verbosity:             int
32         """
33         self.QUERY                     = 0
34         self.REFERENCE                 = 1
35         self.WORKING                   = 2
36         self.INPUTTYPES                = (self.QUERY, self.REFERENCE)
37         self.INPUTWORKINGTYPES         = (self.QUERY, self.REFERENCE, self.WORKING)
38         self.typeToString              = {self.QUERY: "Query", self.REFERENCE: "Reference", self.WORKING: "Working"}
39 
40         self.logHandle                 = logHandle
41         self.verbosity                 = verbosity
42         self.mySqlConnection           = MySqlConnection(self.verbosity-1)
43         self.inputTranscriptContainers = [None, None]
44         self.tableNames                = ["tmpQueryTable_%d" % (random.randint(0, 100000)), "tmpReferenceTable_%d" % (random.randint(0, 100000)), "tmpOutputTable_%d" % (random.randint(0, 100000)), "tmpWorkingTable_%d" % (random.randint(0, 100000))]
45         self.mySqlTranscriptWriters    = [MySqlTranscriptWriter(self.mySqlConnection, name, verbosity-1) for name in self.tableNames]
46         self.writer                    = None
47         self.introns                   = [False, False]
48         self.starts                    = [None, None]
49         self.ends                      = [None, None]
50         self.fivePrimes                = [None, None]
51         self.threePrimes               = [None, None]
52         self.minDistance               = None
53         self.maxDistance               = 0
54         self.minOverlap                = 1
55         self.pcOverlap                 = None
56         self.colinear                  = False
57         self.antisense                 = False
58         self.downstreams               = [False, False]
59         self.upstreams                 = [False, False]
60         self.outputDistance            = False
61         self.absolute                  = False
62         self.strandedDistance          = False
63         self.nbTranscripts             = [None, None]
64         self.nbNucleotides             = [None, None]
65         self.normalization             = False
66         self.included                  = False
67         self.including                 = False
68         self.invert                    = False
69         self.notOverlapping            = False
70         self.splitDifference           = False
71         self.multiple                  = False
72         self.odds                      = False
73         self.overlapResults            = None
74         self.oddResults                = None
75         self.outputContainer           = None
76         self.transcriptsToBeStored     = dict([(type, TranscriptList()) for type in self.INPUTWORKINGTYPES])
77         self.nbPrinted                 = 0
78 
79         self.mySqlConnection.createDatabase()
80 
81         
82     def __del__(self):
83         """
84         Destructor
85         Remove all temporary tables
86         """
87         for type in self.INPUTWORKINGTYPES:
88             self.mySqlTranscriptWriters[type].removeTables()
89         self.mySqlConnection.deleteDatabase()
90         
91     def acceptIntrons(self, type, bool):
92         """
93         Compare transcripts or exons only
94         @param type: whether use query/reference data
95         @type  type: int
96         @param bool: include introns or not
97         @type  bool: boolean
98         """
99         self.introns[type] = bool
00 
01         
02     def restrictToStart(self, type, size):
03         """
04         Restrict a list of transcripts to first nucleotides
05         @param type: whether use query/reference data
06         @type  type: int
07         @param size: the size of the transcript to be considered
08         @type  size: int
09         """
10         self.starts[type]  = size
11         self.introns[type] = False
12         
13         
14     def restrictToEnd(self, type, size):
15         """
16         Restrict a list of transcripts to first nucleotides
17         @param type: whether use query/reference data
18         @type  type: int
19         @param size: the size of the transcript to be considered
20         @type  size: int
21         """
22         self.ends[type]    = size
23         self.introns[type] = False
24         
25         
26     def extendFivePrime(self, type, size):
27         """
28         Extend a list of transcripts by their 5' end
29         @param type: whether use query/reference data
30         @type  type: int
31         @param size: size of the extension
32         @type  size: int
33         """
34         self.fivePrimes[type] = size
35         
36 
37     def extendThreePrime(self, type, size):
38         """
39         Extend the list of query transcripts by their 3' end
40         @param type: whether use query/reference data
41         @type  type: int
42         @param size: size of the extension
43         @type  size: int
44         """
45         self.threePrimes[type] = size
46         
47         
48     def setMinDistance(self, distance):
49         """
50         Set the min distance between two transcripts
51         @param distance: distance
52         @type  distance: int
53         """
54         self.minDistance = distance
55 
56 
57     def setMaxDistance(self, distance):
58         """
59         Set the max distance between two transcripts
60         @param distance: distance
61         @type  distance: int
62         """
63         self.maxDistance = distance
64         
65 
66     def setMinOverlap(self, overlap):
67         """
68         Set the minimum number of nucleotides to declare an overlap
69         @param overlap: minimum number of nucleotides
70         @type  overlap: int
71         """
72         self.minOverlap = overlap
73         
74 
75     def setPcOverlap(self, overlap):
76         """
77         Set the percentage of nucleotides to declare an overlap
78         @param overlap: percentage of nucleotides
79         @type  overlap: int
80         """
81         self.pcOverlap = overlap
82         
83 
84     def setUpstream(self, type, boolean):
85         """
86         Consider transcripts which are upstream of some transcripts
87         @param type:        whether use query/reference data
88         @type  type:        int
89         @param boolean: consider only these transcripts or not
90         @type  boolean: boolean
91         """
92         self.upstreams[type] = boolean
93 
94 
95     def setDownstream(self, type, boolean):
96         """
97         Consider transcripts which are downstream of some transcripts
98         @param type:        whether use query/reference data
99         @type  type:        int
00         @param boolean: consider only these transcripts or not
01         @type  boolean: boolean
02         """
03         self.downstreams[type] = boolean
04 
05 
06     def setOutputDistance(self, boolean):
07         """
08         Output distance between query and reference instead of query transcript
09         @param boolean: whether distance should be output
10         @type  boolean: boolean
11         """
12         self.outputDistance = boolean
13         
14 
15     def setAbsolute(self, boolean):
16         """
17         Do not consider strand when computing distance (thus, having only non-negative values)
18         @param boolean: whether we should consider strands
19         @type  boolean: boolean
20         """
21         self.absolute = boolean
22         
23 
24     def setStrandedDistance(self, boolean):
25         """
26         Return two distance distributions, one per strand
27         @param boolean: whether we should return 2 distance distance
28         @type  boolean: boolean
29         """
30         self.strandedDistance = boolean
31         
32 
33     def getColinearOnly(self, boolean):
34         """
35         Only consider transcripts that overlap in the same direction
36         @param boolean: whether transcripts should overlap in the same direction
37         @type  boolean: boolean
38         """
39         self.colinear = boolean
40         
41                 
42     def getAntisenseOnly(self, boolean):
43         """
44         Only consider transcripts that overlap in the opposite direction
45         @param boolean: whether transcripts should overlap in the opposite direction
46         @type  boolean: boolean
47         """
48         self.antisense = boolean
49         
50         
51     def setIncludedOnly(self, boolean):
52         """
53         Keep the elements from first set which are included in the second set
54         @param boolean: whether to keep included elements only
55         @type  boolean: boolean
56         """
57         self.included = boolean
58 
59 
60     def setIncludingOnly(self, boolean):
61         """
62         Keep the elements from second set which are included in the first set
63         @param boolean: whether to keep included elements only
64         @type  boolean: boolean
65         """
66         self.including = boolean
67 
68 
69     def setNormalization(self, boolean):
70         """
71         Normalize the elements by the number of mappings in the genome
72         @param boolean: whether normalize
73         @type  boolean: boolean
74         """
75         self.normalization = boolean
76 
77 
78     def getInvert(self, boolean):
79         """
80         Only consider transcripts that do not overlap
81         @param boolean: whether invert the selection
82         @type  boolean: boolean
83         """
84         self.invert = boolean
85         
86         
87     def includeNotOverlapping(self, boolean):
88         """
89         Also output the elements which do not overlap
90         @param boolean: whether output the elements which do not overlap
91         @type  boolean: boolean
92         """
93         self.notOverlapping = boolean
94         
95         
96     def setSplitDifference(self, boolean):
97         """
98         Split into intervals when computing difference
99         @param boolean: whether to split
00         @type  boolean: boolean
01         """
02         self.splitDifference = boolean
03 
04 
05     def aggregate(self, boolean):
06         """
07         In merge mode, aggregate multiple transcripts
08         @param boolean: aggregate multiple transcripts
09         @type  boolean: boolean
10         """
11         self.multiple = boolean
12     
13 
14     def getTables(self, type):
15         """
16         Get the SQL tables
17         @param type: type of the table (query, reference, etc.)
18         @type  type: int
19         """
20         return self.mySqlTranscriptWriters[type].getTables()
21     
22 
23     def computeOdds(self, boolean):
24         """
25         Compute odds
26         @param boolean: whether odds should be computed
27         @type  boolean: boolean
28         """
29         self.odds = boolean
30         if self.odds:
31             self.overlapResults = dict()
32         
33         
34     def computeOddsPerTranscript(self, boolean):
35         """
36         Compute odds for each transcript
37         @param boolean: whether odds for each transcript should be computed
38         @type  boolean: boolean
39         """
40         self.odds = boolean
41         if self.odds:
42             self.overlapResults = dict()
43         
44         
45     def removeTables(self):
46         """
47         Remove the temporary MySQL tables
48         """
49         for type in self.INPUTWORKINGTYPES:
50             for chromosome in self.getTables(type):
51                 self.getTables(type)[chromosome].remove()
52 
53 
54     def clearTables(self):
55         """
56         Empty the content of the databases
57         """
58         for type in self.INPUTWORKINGTYPES:
59             if self.transcriptListParsers[type] != None:
60                 for chromosome in self.getTables(type):
61                     self.getTables(type)[chromosome].clear()
62 
63 
64     def extendTranscript(self, type, transcript):
65         """
66         Extend a transcript corresponding to the parameters of the class
67         @param transcript: a transcript
68         @type  transcript: class L{Transcript<Transcript>}
69         @return: the possibly extended transcript
70         """
71         extendedTranscript = Transcript()
72         extendedTranscript.copy(transcript)
73         if self.starts[type] != None:
74             extendedTranscript.restrictStart(self.starts[type])
75         if self.ends[type] != None:
76             extendedTranscript.restrictEnd(self.ends[type])
77         if self.fivePrimes[type] != None:
78             extendedTranscript.extendStart(self.fivePrimes[type])
79         if self.threePrimes[type] != None:
80             extendedTranscript.extendEnd(self.threePrimes[type])
81         return extendedTranscript
82 
83 
84     def storeTranscript(self, type, transcript, now = True):
85         """
86         Add a transcript to a MySQL database, or postpone the store
87         @param type:       whether use query/reference table
88         @type  type:       int
89         @param transcript: a transcript
90         @type  transcript: class L{Transcript<Transcript>}
91         @param now:        whether transcript should be stored now (or stored can be postponed)
92         @type  now:        bool
93         """
94         self.mySqlTranscriptWriters[type].addTranscript(transcript)
95         if type == self.REFERENCE:
96             self.mySqlTranscriptWriters[self.WORKING].addTranscript(transcript)
97         if now:
98             self.mySqlTranscriptWriters[type].write()
99             if type == self.REFERENCE:
00                 self.mySqlTranscriptWriters[self.WORKING].write()
01 
02 
03     def writeTranscript(self, transcript):
04         """
05         Write a transcript in the output file
06         @param transcript: a transcript
07         @type  transcript: class L{Transcript<Transcript>}
08         """
09         if self.writer != None:
10             self.writer.addTranscript(transcript)
11             self.nbPrinted += 1
12 
13 
14     def flushData(self, type = None):
15         """
16         Store the remaining transcripts
17         @param type: whether use query/reference table (None for all)
18         @type  type: int or None
19         """
20         if type == None:
21             types = self.INPUTWORKINGTYPES
22         else:
23             types = [type]
24         for type in types:
25             self.mySqlTranscriptWriters[type].write()
26         if self.writer != None:
27             self.writer.write()
28         
29         
30     def unstoreTranscript(self, type, transcript):
31         """
32         Remove a transcript from a MySQL database
33         @param type:       whether use query/reference table
34         @type  type:       int
35         @param transcript: a transcript
36         @type  transcript: class L{Transcript<Transcript>}
37         """
38         self.getTables(type)[transcript.getChromosome()].removeTranscript(transcript)
39         if type == self.REFERENCE:
40             self.getTables(self.WORKING)[transcript.getChromosome()].removeTranscript(transcript)
41 
42 
43     def addIndexes(self, tables):
44         """
45         Add useful indexes to the tables
46         @param tables: which tables should be indexed
47         @type  tables: list of int
48         """
49         for type in tables:
50             for chromosome in self.getTables(type):
51                 self.getTables(type)[chromosome].createIndex("iStart_transcript_%s_%d_%d" % (chromosome, type, random.randint(0, 100000)), ["start"])
52                 self.getTables(type)[chromosome].exonsTable.createIndex("iTranscriptId_exon_%s_%d_%d" % (chromosome, type, random.randint(0, 100000)), ["transcriptId"])
53 
54 
55     def storeTranscriptList(self, type, transcriptListParser, extension):
56         """
57         Store a transcript list into database
58         @param type:      whether use query/reference parser
59         @type  type:      int
60         @param parser:    a parser of transcript list
61         @type  parser:    class L{TranscriptContainer<TranscriptContainer>}
62         @param extension: extend (or not) the transcripts
63         @type  extension: boolean
64         """
65         progress = Progress(transcriptListParser.getNbTranscripts(), "Writing transcripts for %s" % ("query" if type == self.QUERY else "reference"), self.verbosity-1)
66         for transcript in transcriptListParser.getIterator():
67             if extension:
68                 transcript = self.extendTranscript(type, transcript)
69             self.mySqlTranscriptWriters[type].addTranscript(transcript)
70             progress.inc()
71         self.mySqlTranscriptWriters[type].write()
72         progress.done()
73         if type == self.REFERENCE:
74             for chromosome in self.getTables(self.REFERENCE):
75                 self.getTables(self.WORKING)[chromosome] = MySqlTranscriptTable(self.mySqlConnection, self.tableNames[self.WORKING], chromosome, self.verbosity-1)
76                 self.getTables(self.WORKING)[chromosome].copy(self.getTables(self.REFERENCE)[chromosome])
77         
78     
79     def setInputTranscriptContainer(self, type, inputTranscriptContainer):
80         """
81         Set an input transcript list container
82         @param type:                     whether use query/reference parser
83         @type  type:                     int
84         @param inputTranscriptContainer: a container
85         @type  inputTranscriptContainer: class L{TranscriptContainer<TranscriptContainer>}
86         """
87         self.inputTranscriptContainers[type] = inputTranscriptContainer
88         self.nbTranscripts[type]             = self.inputTranscriptContainers[type].getNbTranscripts()
89         self.nbNucleotides[type]             = self.inputTranscriptContainers[type].getNbNucleotides()
90 
91 
92     def setOutputWriter(self, writer):
93         """
94         Set an output transcript list writer
95         @param writer: a writer
96         @type  writer: class L{TranscriptListWriter<TranscriptListWriter>}
97         """
98         self.writer = writer
99 
00 
01     def compareTranscript(self, transcript1, transcript2, includeDistance = False):
02         """
03         Compare two transcripts, using user defined parameters
04         @param transcript1:     a transcript from the query set (already extended)
05         @type  transcript1:     class L{Transcript<Transcript>}
06         @param transcript2:     a transcript from the reference set (already extended)
07         @type  transcript2:     class L{Transcript<Transcript>}
08         @param includeDistance: take into account the distance too
09         @type  includeDistance: boolean
10         @return:                true, if they overlap
11         """
12         extendedTranscript1 = Transcript()
13         extendedTranscript1.copy(transcript1)
14         if includeDistance:
15             if self.maxDistance > 0:
16                 extendedTranscript1.extendStart(self.maxDistance)
17                 extendedTranscript1.extendEnd(self.maxDistance)
18 
19         minOverlap = self.minOverlap
20         if self.pcOverlap != None:
21             minOverlap = max(minOverlap, transcript1.getSize() / 100.0 * self.pcOverlap)
22         if not extendedTranscript1.overlapWith(transcript2, self.minOverlap):
23             return False
24         if (self.downstreams[self.QUERY]       and transcript2.getStart() > extendedTranscript1.getStart()) or \
25              (self.upstreams[self.QUERY]       and transcript2.getEnd()   < extendedTranscript1.getEnd())   or \
26              (self.downstreams[self.REFERENCE] and extendedTranscript1.getStart() > transcript2.getStart()) or \
27              (self.upstreams[self.REFERENCE]   and extendedTranscript1.getEnd()   < transcript2.getEnd()):
28             return False
29         if (self.antisense and extendedTranscript1.getDirection() == transcript2.getDirection()) or (self.colinear and extendedTranscript1.getDirection() != transcript2.getDirection()):
30             return False
31         if self.included and not transcript2.include(extendedTranscript1):
32             return False
33         if self.including and not extendedTranscript1.include(transcript2):
34             return False
35         if self.introns[self.REFERENCE] and self.introns[self.QUERY]:
36             if self.logHandle != None:
37                 self.logHandle.write("%s overlaps with intron of %s\n" % (str(extendedTranscript1), str(transcript2)))
38             return True
39         if (not self.introns[self.REFERENCE]) and (not self.introns[self.QUERY]) and extendedTranscript1.overlapWithExon(transcript2, minOverlap):
40             if self.logHandle != None:
41                 self.logHandle.write("%s overlaps with exon of %s\n" % (str(extendedTranscript1), str(transcript2)))
42             return True
43         return False
44 
45 
46     def compareTranscriptToList(self, transcript1):
47         """
48         Compare a transcript to the reference list of transcripts
49         (Do not extend the transcripts, except for the distance)
50         @param transcript1: a transcript (from the query set)
51         @type  transcript1: class L{Transcript<Transcript>}
52         @return:            the reference transcripts overlapping
53         """
54         # no transcript in the reference table
55         if transcript1.getChromosome() not in self.getTables(self.WORKING):
56             return
57 
58         # retrieve the the transcripts that may overlap in the working tables
59         clauses = []
60         extendedTranscript1 = Transcript()
61         extendedTranscript1.copy(transcript1)
62         if self.maxDistance > 0:
63             extendedTranscript1.extendStart(self.maxDistance)
64         if self.maxDistance > 0:
65             extendedTranscript1.extendEnd(self.maxDistance)
66         command = "SELECT * FROM %s WHERE (" % (self.getTables(self.WORKING)[transcript1.getChromosome()].getName())
67         for binPair in extendedTranscript1.getBins():
68             clause = "bin "
69             if binPair[0] == binPair[1]:
70                 clause += "= %i" % (binPair[0])
71             else:
72                 clause += "BETWEEN %i AND %i" % (binPair[0], binPair[1])
73             clauses.append(clause)
74         command += " OR ".join(clauses)
75         command += ") AND start <= %d AND end >= %d" % (extendedTranscript1.getEnd(), extendedTranscript1.getStart())
76 
77         for index2, transcript2 in self.getTables(self.REFERENCE)[transcript1.getChromosome()].selectTranscripts(command):
78             if self.compareTranscript(extendedTranscript1, transcript2):
79                 yield transcript2
80 
81     
82     def compareTranscriptList(self):
83         """
84         Compare a list of transcript to the reference one
85         @return: the transcripts that overlap with the reference set
86         """
87         distance      = 0
88         nbClustersIn  = 0
89         nbClustersOut = 0
90         if self.maxDistance != None:
91             distance = self.maxDistance
92 
93         self.addIndexes([self.QUERY, self.REFERENCE])
94 
95         # export the container into tables
96         self.storeTranscriptList(self.QUERY, self.inputTranscriptContainers[self.QUERY], True)
97         self.storeTranscriptList(self.REFERENCE, self.inputTranscriptContainers[self.REFERENCE], True)
98 
99         # looping
00         for chromosome1 in sorted(self.getTables(self.QUERY).keys()):
01             # get range of transcripts
02             command = "SELECT MIN(start), MAX(end), COUNT(id) FROM %s" % (self.getTables(self.QUERY)[chromosome1].getName())
03             query     = self.mySqlConnection.executeQuery(command)
04             result    = query.getLine()
05             first     = result[0]
06             last      = result[1]
07             nb        = result[2]
08 
09             transcripts1 = []
10             toBeRemoved1 = []
11             transcripts2 = []
12             toBeRemoved2 = []
13             overlapsWith = []
14             nbOverlaps   = []
15             nbChunks     = max(1, nb / 100)
16             chunkSize    = (last - first) / nbChunks
17             progress     = Progress(nbChunks + 1, "Analyzing chromosome %s" % (chromosome1), self.verbosity-1)
18             for chunk in range(nbChunks + 1):
19 
20                 # load transcripts
21                 start   = first + chunk * chunkSize
22                 end     = start + chunkSize - 1
23                 command = "SELECT * FROM %s WHERE start BETWEEN %d AND %d" % (self.getTables(self.QUERY)[chromosome1].getName(), start, end-1)
24                 for index1, transcript1 in self.getTables(self.QUERY)[chromosome1].selectTranscripts(command):
25                     transcripts1.append(transcript1)
26                     overlapsWith.append([])
27                     nbOverlaps.append(0)
28                     nbClustersIn += 1 if "nbElements" not in transcript1.getTagNames() else transcript1.getTagValue("nbElements")
29                 command = "DELETE FROM %s WHERE start < %d" % (self.getTables(self.QUERY)[chromosome1].getName(), end)
30                 self.mySqlConnection.executeQuery(command)
31                 if chromosome1 in self.getTables(self.REFERENCE):
32                     command = "SELECT * FROM %s WHERE start BETWEEN %d AND %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), start-distance, end+distance-1)
33                     if chunk == 0:
34                         command = "SELECT * FROM %s WHERE start < %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), end + distance)
35                     for index2, transcript2 in self.getTables(self.REFERENCE)[chromosome1].selectTranscripts(command):
36                         transcripts2.append(transcript2)
37                     command = "DELETE FROM %s WHERE start < %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), end + distance)
38                     self.mySqlConnection.executeQuery(command)
39 
40                 # compare sets
41                 for index1, transcript1 in enumerate(transcripts1):
42                     overlappingNames = []
43                     nbElements1      = 1 if "nbElements" not in transcript1.getTagNames() else transcript1.getTagValue("nbElements")
44                     for transcript2 in transcripts2:
45                         if self.compareTranscript(transcript1, transcript2, True):
46                             id2 = transcript2.getTagValue("ID") if "ID" in transcript2.getTagNames() else transcript2.getName()
47                             if id2 not in overlapsWith[index1]:
48                                 overlapsWith[index1].append(id2)
49                                 nbOverlaps[index1] += 1 if "nbElements" not in transcript2.getTagNames() else transcript2.getTagValue("nbElements")
50                                 if self.odds:
51                                     if transcript2.getName() not in self.overlapResults:
52                                         self.overlapResults[transcript2.getName()] = 0
53                                     self.overlapResults[transcript2.getName()] += nbElements1
54 
55                     # check if query transcript extends bounds of the chunk
56                     if transcript1.getEnd() < end:
57                         if Utils.xor(overlapsWith[index1], self.invert) or self.notOverlapping:
58                             if overlapsWith[index1]:
59                                 transcript1.setTagValue("overlapWith", ",".join(overlapsWith[index1])[:100])
60                                 transcript1.setTagValue("nbOverlaps", "%d" % (nbOverlaps[index1]))
61                             elif self.notOverlapping:
62                                 transcript1.setTagValue("nbOverlaps", "0")
63                             self.writeTranscript(transcript1)
64                             nbClustersOut += nbElements1
65                         toBeRemoved1.append(index1)
66 
67                 # update list of query transcripts
68                 for index1 in reversed(toBeRemoved1):
69                     del transcripts1[index1]
70                     del overlapsWith[index1]
71                     del nbOverlaps[index1]
72                 toBeRemoved1 = []
73 
74                 # check if the reference transcripts extends bounds of the chunk
75                 for index2, transcript2 in enumerate(transcripts2):
76                     if transcript2.getEnd() + distance < end:
77                         toBeRemoved2.append(index2)
78                 for index2 in reversed(toBeRemoved2):
79                     del transcripts2[index2]
80                 toBeRemoved2 = []
81 
82                 progress.inc()
83 
84             for index1, transcript1 in enumerate(transcripts1):
85                 if Utils.xor(overlapsWith[index1], self.invert) or self.notOverlapping:
86                     if overlapsWith[index1]:
87                         transcript1.setTagValue("overlapWith", ",".join(overlapsWith[index1])[:100])
88                         transcript1.setTagValue("nbOverlaps", "%d" % (nbOverlaps[index1]))
89                     elif self.notOverlapping:
90                         transcript1.setTagValue("nbOverlaps", "0")
91                     self.writeTranscript(transcript1)
92                     nbClustersOut += 1 if "nbElements" not in transcript1.getTagNames() else transcript1.getTagValue("nbElements")
93             progress.done()
94             self.getTables(self.QUERY)[chromosome1].remove()
95             if chromosome1 in self.getTables(self.REFERENCE):
96                 self.getTables(self.REFERENCE)[chromosome1].remove()
97                 self.getTables(self.WORKING)[chromosome1].remove()
98 
99         self.flushData()
00         if self.writer != None:
01             self.writer.close()
02             self.writer = None
03 
04         if self.verbosity > 0:
05             print "reference: %d elements" % (self.nbTranscripts[self.REFERENCE])
06             print "query:     %d elements, %d clustered" % (self.nbTranscripts[self.QUERY], nbClustersIn)
07             if self.nbTranscripts[self.QUERY] != 0:
08                 print "output:    %d elements (%.2f%%)"% (self.nbPrinted, self.nbPrinted / float(self.nbTranscripts[self.QUERY]) * 100), 
09                 if nbClustersOut != 0:
10                     print ", %d clustered (%.2f%%)" % (nbClustersOut, float(nbClustersOut) / nbClustersIn * 100)
11 
12 
13     def compareTranscriptListDistance(self):
14         """
15         Compare a list of transcript to the reference one
16         @return: the distance distributions in a hash
17         """
18         nbDistances  = 0
19         distances    = {}
20         absDistances = {}
21         strandedDistances = dict([(strand, {}) for strand in (1, -1)])
22         
23         # export the container into tables
24         self.storeTranscriptList(self.QUERY, self.inputTranscriptContainers[self.QUERY], True)
25         self.storeTranscriptList(self.REFERENCE, self.inputTranscriptContainers[self.REFERENCE], True)
26         
27         progress = Progress(self.nbTranscripts[self.QUERY], "Analyzing chromosomes", self.verbosity-1)
28         for transcript1 in self.inputTranscriptContainers[self.QUERY].getIterator():
29             # get the distance
30             transcript1    = self.extendTranscript(self.QUERY, transcript1)
31             distance       = self.maxDistance + 1
32             strand         = None
33             closestElement = "None"
34             for transcript2 in self.compareTranscriptToList(transcript1):
35                 thisStrand = transcript1.getDirection() * transcript2.getDirection()
36                 if self.antisense or (not self.colinear and transcript1.getDirection() != transcript2.getDirection()):
37                     transcript2.reverse()
38                 if self.absolute:
39                     transcript2.setDirection(transcript1.getDirection())
40                 if transcript2.getDirection() == transcript1.getDirection():
41                     if self.starts[self.REFERENCE] != None:
42                         transcript2.restrictStart(self.starts[self.REFERENCE])
43                     if self.ends[self.REFERENCE] != None:
44                         transcript2.restrictEnd(self.ends[self.REFERENCE])
45                     thisDistance = transcript1.getRelativeDistance(transcript2)
46                     if (self.absolute):
47                         thisDistance = abs(thisDistance)
48                     if abs(thisDistance) < abs(distance):
49                         distance       = thisDistance
50                         strand         = thisStrand
51                         closestElement = transcript2.getTagValue("ID") if "ID" in transcript2.getTagNames() else transcript2.getName()
52             if (distance <= self.maxDistance) and (self.minDistance == None or distance >= self.minDistance):
53                 nbDistances                        += 1
54                 distances[distance]                 = distances.get(distance, 0) + 1
55                 absDistance                         = abs(distance)
56                 absDistances[absDistance]           = absDistances.get(absDistance, 0) + 1
57                 strandedDistances[strand][distance] = strandedDistances[strand].get(distance, 0)
58                 if distance not in strandedDistances[-strand]:
59                     strandedDistances[-strand][distance] = 0
60 
61             # write transcript
62             if distance == self.maxDistance + 1:
63                 distance = "None"
64             tmpTranscript = Transcript()
65             tmpTranscript.copy(transcript1)
66             tmpTranscript.setTagValue("distance", distance)
67             tmpTranscript.setTagValue("closestElement", closestElement)
68             self.writeTranscript(tmpTranscript)
69 
70             progress.inc()
71         progress.done()
72 
73         self.flushData()
74 
75         if self.verbosity > 0:
76             print "reference: %d sequences" % (self.nbTranscripts[self.REFERENCE])
77             print "query:     %d sequences" % (self.nbTranscripts[self.QUERY])
78             if nbDistances == 0:
79                 print "Nothing matches"
80             else:
81                 print "min/avg/med/max transcripts: %d/%.2f/%.1f/%d" % Utils.getMinAvgMedMax(absDistances)
82                 print "for %d distances (%.2f%%)" % (nbDistances, float(nbDistances) / self.nbTranscripts[self.QUERY] * 100)
83 
84         if self.strandedDistance:
85             return strandedDistances
86         return distances
87 
88 
89     def compareTranscriptListMerge(self):
90         """
91         Merge the query list of transcript with itself
92         @return: the merged transcripts in a transcript list database
93         """
94         nbMerges  = 0
95 
96         for type in (self.QUERY, self.REFERENCE):
97             self.storeTranscriptList(type, self.inputTranscriptContainers[type], True)
98         self.flushData()
99 
00         # Loop on the chromosomes
01         for chromosome in sorted(self.getTables(self.QUERY).keys()):
02             if chromosome not in self.getTables(self.REFERENCE):
03                 continue
04 
05             # Get the size of the chromosome
06             maxEnd   = 0
07             nbChunks = 0
08             for type in (self.QUERY, self.REFERENCE):
09                 command  = "SELECT MAX(end) from %s" % (self.getTables(type)[chromosome].getName())
10                 query    = self.mySqlConnection.executeQuery(command)
11                 maxEnd   = max(maxEnd, int(query.getLine()[0]))
12                 nbChunks = max(nbChunks, self.getTables(type)[chromosome].getNbElements())
13 
14             mergedTranscripts = {}
15             transcripts = {self.QUERY: [], self.REFERENCE: []}
16             progress    = Progress(nbChunks, "Analyzing %s" % (chromosome), self.verbosity-1)
17             for i in range(nbChunks):
18                 rangeStart = int(i * (float(maxEnd) / nbChunks)) + 1
19                 rangeEnd   = int((i+1) * (float(maxEnd) / nbChunks))
20 
21                 # Get all transcripts in query and reference from chunk
22                 for type in (self.QUERY, self.REFERENCE):
23                     correction = 0 if self.QUERY else self.maxDistance
24                     command = "SELECT * FROM %s WHERE start <= %d" % (self.getTables(type)[chromosome].getName(), rangeEnd + correction)
25                     for index, transcript in self.getTables(type)[chromosome].selectTranscripts(command):
26                         transcripts[type].append(transcript)
27 
28                 # Merge elements between the two samples
29                 for iQuery, queryTranscript in enumerate(transcripts[self.QUERY]):
30                     for iReference, referenceTranscript in enumerate(transcripts[self.REFERENCE]):
31                         if referenceTranscript == None: continue
32                         if self.compareTranscript(queryTranscript, referenceTranscript, True):
33                             if queryTranscript.getDirection() != referenceTranscript.getDirection():
34                                 referenceTranscript.setDirection(queryTranscript.getDirection())
35                             queryTranscript.merge(referenceTranscript, self.normalization)
36                             nbMerges += 1
37                             transcripts[self.REFERENCE][iReference] = None
38                             if not self.multiple:
39                                 mergedTranscripts[iQuery] = 0
40 
41                 # Remove transcripts from database
42                 for type in (self.QUERY, self.REFERENCE):
43                     correction = 0 if self.QUERY else self.maxDistance
44                     command    = "DELETE FROM %s WHERE start <= %d" % (self.getTables(type)[chromosome].getName(), rangeEnd - correction)
45                     query      = self.mySqlConnection.executeQuery(command)
46 
47                 # Just in case, self-merge the elements in the query (beware of mergedTranscripts!)
48                 if (self.multiple):
49                     for iQuery1, queryTranscript1 in enumerate(transcripts[self.QUERY]):
50                         if queryTranscript1 == None: continue
51                         for iQuery2, queryTranscript2 in enumerate(transcripts[self.QUERY]):
52                             if iQuery2 <= iQuery1 or queryTranscript2 == None: continue
53                             minOverlap = self.minOverlap
54                             if self.pcOverlap != None:
55                                 minOverlap = max(minOverlap, queryTranscript1.getSize() / 100.0 * self.pcOverlap)
56                             if queryTranscript2.overlapWith(queryTranscript1, minOverlap) and (queryTranscript1.getDirection() == queryTranscript2.getDirection() or not self.colinear):
57                                 if queryTranscript1.getDirection() != queryTranscript2.getDirection():
58                                     queryTranscript2.setDirection(queryTranscript1.getDirection())
59                                 queryTranscript1.merge(queryTranscript2, self.normalization)
60                                 transcripts[self.QUERY][iQuery2] = None
61                                 nbMerges += 1
62                                 if not self.multiple:
63                                     mergedTranscripts[iQuery1] = 0
64 
65                 # Update the sets of transcripts and write into database (also update mergedTranscripts)
66                 newTranscripts = {self.QUERY: [], self.REFERENCE: []}
67                 newMergedTranscripts = {}
68                 for type in (self.QUERY, self.REFERENCE):
69                     for i, transcript in enumerate(transcripts[type]):
70                         if transcript == None: continue
71                         correction = 0 if self.QUERY else self.maxDistance
72                         if transcript.getEnd() < rangeEnd - correction:
73                             if self.multiple or ((type == self.QUERY) and (i in mergedTranscripts)):
74                                 self.writeTranscript(transcripts[type][i])
75                         else:
76                             if type == self.QUERY and i in mergedTranscripts:
77                                 newMergedTranscripts[len(newTranscripts[type])] = 0
78                             newTranscripts[type].append(transcript)
79                 transcripts = newTranscripts
80                 mergedTranscripts = newMergedTranscripts
81 
82                 progress.inc()
83             progress.done()
84 
85             for type in (self.QUERY, self.REFERENCE):
86                 for i, transcript in enumerate(transcripts[type]):
87                     if transcripts == None: continue
88                     if self.multiple or ((type == self.QUERY) and (i in mergedTranscripts)):
89                         self.writeTranscript(transcripts[type][i])
90 
91         # Manage chromosomes with no corresponding data
92         if self.multiple:
93             for type in self.INPUTTYPES:
94                 for chromosome in self.getTables(type):
95                     if chromosome in self.getTables(1 - type):
96                         continue
97                     for transcript in self.getTables(self.OUTPUT)[chromosome].getIterator():
98                         self.writeTranscript(transcript)
99 
00         self.flushData()
01         if self.writer != None:
02             self.writer.close()
03             self.writer = None
04 
05         if self.verbosity > 0:
06             print "query:    %d sequences" % (self.nbTranscripts[self.QUERY])
07             print "# merges: %d" % (nbMerges)
08             print "# printed %d (%.2f%%)" % (self.nbPrinted, self.nbPrinted / float(self.nbTranscripts[self.QUERY]) * 100)
09 
10 
11     def compareTranscriptListSelfMerge(self):
12         """
13         Merge the query list of transcript with itself
14         @return: the merged transcripts in a transcript list database
15         """
16         nbMerges  = 0
17         distance  = self.maxDistance if self.maxDistance != None else 0
18 
19         self.addIndexes([self.QUERY])
20         self.storeTranscriptList(self.QUERY, self.inputTranscriptContainers[self.QUERY], True)
21         self.flushData()
22 
23         # looping
24         for chromosome1 in sorted(self.getTables(self.QUERY).keys()):
25             transcripts2 = []
26 
27             # get range of transcripts
28             progress = Progress(self.getTables(self.QUERY)[chromosome1].getNbElements(), "Analyzing chromosome %s" % (chromosome1), self.verbosity-1)
29             command  = "SELECT * FROM %s ORDER BY start" % (self.getTables(self.QUERY)[chromosome1].getName())
30             for index1, transcript1 in self.getTables(self.QUERY)[chromosome1].selectTranscripts(command):
31 
32                 # compare sets
33                 toBeRemoved = set()
34                 toBePrinted = set()
35                 for index2, transcript2 in enumerate(transcripts2):
36 
37                     if self.compareTranscript(transcript1, transcript2, True):
38                         if transcript1.getDirection() != transcript2.getDirection():
39                             transcript2.setDirection(transcript1.getDirection())
40                         transcript1.merge(transcript2, self.normalization)
41                         toBeRemoved.add(index2)
42                         nbMerges += 1
43                     elif transcript2.getEnd() + distance < transcript1.getStart():
44                         toBePrinted.add(index2)
45                 transcripts2.append(transcript1)
46 
47                 for index2 in sorted(toBePrinted):
48                     self.writeTranscript(transcripts2[index2])
49                 transcripts2 = [transcripts2[index2] for index2 in range(len(transcripts2)) if index2 not in (toBeRemoved | toBePrinted)]
50 
51             for transcript2 in transcripts2:
52                 self.writeTranscript(transcript2)
53             progress.done()
54             self.getTables(self.QUERY)[chromosome1].remove()
55 
56         self.flushData()
57         if self.writer != None:
58             self.writer.close()
59             self.writer = None
60 
61         if self.verbosity > 0:
62             print "query:    %d sequences" % (self.nbTranscripts[self.QUERY])
63             print "# merges: %d" % (nbMerges)
64             print "# printed %d (%.2f%%)" % (self.nbPrinted, self.nbPrinted / float(self.nbTranscripts[self.QUERY]) * 100)
65 
66 
67     def getDifferenceTranscriptList(self):
68         """
69         Get the elements of the first list which do not overlap the second list (at the nucleotide level)
70         @return: the transcripts that overlap with the reference set
71         """
72         distance = 0 if self.maxDistance == None else self.maxDistance
73 
74         self.addIndexes([self.QUERY, self.REFERENCE])
75 
76         # export the container into tables
77         self.storeTranscriptList(self.QUERY, self.inputTranscriptContainers[self.QUERY], True)
78         self.storeTranscriptList(self.REFERENCE, self.inputTranscriptContainers[self.REFERENCE], True)
79 
80         # looping
81         for chromosome1 in sorted(self.getTables(self.QUERY).keys()):
82             # get range of transcripts
83             command = "SELECT MIN(start), MAX(end), COUNT(id) FROM %s" % (self.getTables(self.QUERY)[chromosome1].getName())
84             query   = self.mySqlConnection.executeQuery(command)
85             result  = query.getLine()
86             first   = result[0]
87             last    = result[1]
88             nb      = result[2]
89 
90             transcripts1 = []
91             transcripts2 = []
92             nbChunks     = max(1, nb / 100)
93             chunkSize    = (last - first) / nbChunks
94             progress     = Progress(nbChunks + 1, "Analyzing chromosome %s" % (chromosome1), self.verbosity-1)
95             for chunk in range(nbChunks + 1):
96 
97                 # load transcripts
98                 start   = first + chunk * chunkSize
99                 end     = start + chunkSize - 1
00                 command = "SELECT * FROM %s WHERE start BETWEEN %d AND %d" % (self.getTables(self.QUERY)[chromosome1].getName(), start, end-1)
01                 for index1, transcript1 in self.getTables(self.QUERY)[chromosome1].selectTranscripts(command):
02                     transcripts1.append(transcript1)
03                 command = "DELETE FROM %s WHERE start < %d" % (self.getTables(self.QUERY)[chromosome1].getName(), end)
04                 self.mySqlConnection.executeQuery(command)
05                 if chromosome1 in self.getTables(self.REFERENCE):
06                     command = "SELECT * FROM %s WHERE start BETWEEN %d AND %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), start-distance, end+distance-1)
07                     if chunk == 0:
08                         command = "SELECT * FROM %s WHERE start < %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), end + distance)
09                     for index2, transcript2 in self.getTables(self.REFERENCE)[chromosome1].selectTranscripts(command):
10                         transcripts2.append(transcript2)
11                     command = "DELETE FROM %s WHERE start < %d" % (self.getTables(self.REFERENCE)[chromosome1].getName(), end + distance)
12                     self.mySqlConnection.executeQuery(command)
13 
14                 # compare sets
15                 toBeRemoved1 = []
16                 for index1, transcript1 in enumerate(transcripts1):
17                     newTranscript1 = Transcript()
18                     newTranscript1.copy(transcript1)
19                     for transcript2 in transcripts2:
20                         newTranscript1 = newTranscript1.getDifference(transcript2)
21                         if newTranscript1 == None:
22                             toBeRemoved1.append(index1)
23                             break
24                     transcripts1[index1] = newTranscript1
25 
26                     # check if query transcript extends bounds of the chunk
27                     if newTranscript1 != None and newTranscript1.getEnd() < end:
28                         if self.splitDifference:
29                             for exon in newTranscript1.getExons():
30                                 transcript = Transcript()
31                                 transcript.copy(exon)
32                                 self.writeTranscript(transcript)
33                         else:
34                             self.writeTranscript(newTranscript1)
35                         toBeRemoved1.append(index1)
36 
37                 # update list of query transcripts
38                 for index1 in reversed(toBeRemoved1):
39                     del transcripts1[index1]
40 
41                 # check if the reference transcripts extends bounds of the chunk
42                 toBeRemoved2 = []
43                 for index2, transcript2 in enumerate(transcripts2):
44                     if transcript2.getEnd() + distance < end:
45                         toBeRemoved2.append(index2)
46                 for index2 in reversed(toBeRemoved2):
47                     del transcripts2[index2]
48 
49                 progress.inc()
50 
51             for transcript1 in transcripts1:
52                 if self.splitDifference:
53                     for exon in transcript1.getExons():
54                         transcript = Transcript()
55                         transcript.copy(exon)
56                         self.writeTranscript(transcript)
57                 else:
58                     self.writeTranscript(transcript1)
59             progress.done()
60             self.getTables(self.QUERY)[chromosome1].remove()
61             if chromosome1 in self.getTables(self.REFERENCE):
62                 self.getTables(self.REFERENCE)[chromosome1].remove()
63                 self.getTables(self.WORKING)[chromosome1].remove()
64 
65         self.flushData()
66         if self.writer != None:
67             self.writer.close()
68             self.writer = None
69 
70         if self.verbosity > 0:
71             print "query:     %d elements" % (self.nbTranscripts[self.QUERY])
72             print "reference: %d elements" % (self.nbTranscripts[self.REFERENCE])
73             print "# printed: %d (%.2f%%)" % (self.nbPrinted, self.nbPrinted / float(self.nbTranscripts[self.QUERY]) * 100)
74 
75 
76     def getOddsPerTranscript(self):
77         """
78         Return overlap results
79         @return a dict of data
80         """
81         if not self.odds:
82             raise Exception("Did not compute odds!")
83         return self.overlapResults
84 
85 
86     def getOdds(self):
87         """
88         Return odds about the overlap
89         @return a dict of data
90         """
91         if not self.odds:
92             raise Exception("Did not compute odds!")
93         if self.oddResults != None:
94             return self.oddResults
95         self.oddResults = {}
96         for name, value in self.overlapResults.iteritems():
97             self.oddResults[value] = self.oddResults.get(value, 0) + 1
98         return self.oddResults