← Back to branch summary

~ubuntu-branches/ubuntu/natty/python-cogent/natty

« back to all changes in this revision

Viewing changes to cogent/draw/dotplot.py

  • Committer: Bazaar Package Importer
  • Author(s): Steffen Moeller
  • Date: 2010-12-04 22:30:35 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20101204223035-j11kinhcrrdgg2p2
Tags: 1.5-1
* Bumped standard to 3.9.1, no changes required.
* New upstream version.
  - major additions to Cookbook
  - added AlleleFreqs attribute to ensembl Variation objects.
  - added getGeneByStableId method to genome objects.
  - added Introns attribute to Transcript objects and an Intron class.
  - added Mann-Whitney test and a Monte-Carlo version
  - exploratory and confirmatory period estimation techniques (suitable for
    symbolic and continuous data)
  - Information theoretic measures (AIC and BIC) added
  - drawing of trees with collapsed nodes
  - progress display indicator support for terminal and GUI apps
  - added parser for illumina HiSeq2000 and GAiix sequence files as 
    cogent.parse.illumina_sequence.MinimalIlluminaSequenceParser.
  - added parser to FASTQ files, one of the output options for illumina's
    workflow, also added cookbook demo.
  - added functionality for parsing of SFF files without the Roche tools in
    cogent.parse.binary_sff
  - thousand fold performance improvement to nmds
  - >10-fold performance improvements to some Table operations

Show diffs side-by-side

added added

removed removed

Lines of Context:
cogent/draw/dotplot.py
1
1
#!/usr/bin/env python
2
2
from __future__ import division
3
 
import logging
4
 
import warnings
5
3
 
6
4
from matplotlib.path import Path
7
5
from matplotlib.patches import PathPatch
8
6
 
 
7
from cogent.util.warning import discontinued
9
8
from cogent.draw.linear import Display
10
9
from cogent.draw.rlg2mpl import Drawable, figureLayout
11
10
from cogent.align.align import dotplot
14
13
__copyright__ = "Copyright 2007-2009, The Cogent Project"
15
14
__credits__ = ["Gavin Huttley", "Peter Maxwell"]
16
15
__license__ = "GPL"
17
 
__version__ = "1.4.1"
 
16
__version__ = "1.5.0"
18
17
__maintainer__ = "Gavin Huttley"
19
18
__email__ = "gavin.huttley@anu.edu.au"
20
19
__status__ = "Production"
21
20
 
22
 
LOG = logging.getLogger("cogent.draw")
23
 
 
24
21
def suitable_threshold(window, desired_probability):
25
22
    """Use cumulative binomial distribution to find the number of identical
26
23
    bases which we expect a nucleotide window-mer to have with the desired
48
45
    
49
46
    Returns a matplotlib axes object placed and scaled ready for plotting 
50
47
    a sequence vs sequence comparison between the sequences (or alignments) 
51
 
    seq1 and seq2, which are also displayed. The longest dimension of the 
52
 
    figure will be inches + 2*margin and its aspect ratio will depend on the 
53
 
    sequence lengths as the sequences are drawn to the same scale"""
 
48
    seq1 and seq2, which are also displayed. The aspect ratio will depend on 
 
49
    the sequence lengths as the sequences are drawn to the same scale"""
54
50
    
55
51
    import matplotlib.pyplot as plt
56
52
    
57
 
    if not isinstance(seq1, Display): seq1 = Display(seq1)
58
 
    if not isinstance(seq2, Display): seq2 = Display(seq2)
59
 
    
60
53
    (x1, y1, w1, h1) = _reinchify(*seq1.figureLayout(
61
54
        labeled=True, bottom=bottom, margin=0))
62
55
    (x2, y2, w2, h2) = _reinchify(*seq2.figureLayout(
86
79
    return d
87
80
 
88
81
class Display2D(Drawable):
89
 
    def __init__(self, seq1, seq2):
90
 
        self.seq1 = seq1
91
 
        self.seq2 = seq2
 
82
    def __init__(self, seq1, seq2, **kw):
 
83
        if not isinstance(seq1, Display):
 
84
            seq1 = Display(seq1, **kw)
 
85
        if not isinstance(seq2, Display):
 
86
            seq2 = Display(seq2, **kw)
 
87
        self.seq1 = seq1.base
 
88
        self.seq1d = seq1
 
89
        self.seq2 = seq2.base
 
90
        self.seq2d = seq2
92
91
        self._cache = {}
 
92
        # Check inputs are sufficiently sequence-like
 
93
        assert len(self.seq1) == len(str(self.seq1))
 
94
        assert len(self.seq2) == len(str(self.seq2))
93
95
    
94
96
    def _calc_lines(self, window, threshold, min_gap):
95
97
        # Cache dotplot line segment coordinates as they can sometimes
99
101
            universe = (len1-window) * (len2-window)
100
102
            acceptable_noise = min(len1, len2) / window
101
103
            threshold = suitable_threshold(window, acceptable_noise/universe)
102
 
            LOG.info('require %s / %s bases' % (threshold, window))
103
 
            LOG.info('expect %s / %s matching' % (acceptable_noise, universe))
 
104
            # print 'require %s / %s bases' % (threshold, window)
 
105
            # print 'expect %s / %s matching' % (acceptable_noise, universe)
104
106
        
105
107
        key = (min_gap, window, threshold)
106
108
        if not self._cache.has_key(key):
107
 
            fwd = dotplot(self.seq1._seq, self.seq2._seq,
108
 
                    window, threshold, min_gap, None, False)
 
109
            fwd = dotplot(str(self.seq1), str(self.seq2),
 
110
                    window, threshold, min_gap, None)
109
111
            if hasattr(self.seq1, "reversecomplement"):
110
 
                rev = dotplot(self.seq1.reversecomplement()._seq, 
111
 
                        self.seq2._seq, window, threshold, min_gap, None, False)
 
112
                rev = dotplot(str(self.seq1.reversecomplement()), 
 
113
                        str(self.seq2), window, threshold, min_gap, None)
112
114
                rev = [((len1-x1,y1),(len1-x2,y2)) for ((x1,y1),(x2,y2)) in rev]
113
115
            else:
114
116
                rev = []
121
123
        # hard to pick min_gap without knowing pixels per base, and
122
124
        # matplotlib is reasonably fast anyway, so:
123
125
        if join_gaps is not None:
124
 
            warnings.warn('"join_gaps" no longer does anything', 
125
 
                    DeprecationWarning, stacklevel=2)
126
 
        ax = comparison_display(self.seq1, self.seq2, **kw)
 
126
            discontinued('argument', 'join_gaps', '1.6')
 
127
        ax = comparison_display(self.seq1d, self.seq2d, **kw)
127
128
        (fwd, rev) = self._calc_lines(window, None, min_gap)
128
 
        LOG.info('lines %s fwd, %s rev' % (len(fwd), len(rev)))
129
129
        for (lines, colour) in [(fwd, 'blue'), (rev, 'red')]:
130
130
            vertices = []
131
131
            for segment in lines: